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(57) Abstract 

Methods are presented 
for calculating the wavelet 
filters and their inverses 
which rely on a new method 
for sampling (UA, TA. A) 
either digital or analog data. 
These methods combine 
and extend to give novel 
procedures for non-reversible 
multi-dimensional data 
compression. For selected 
applications, this procedure 
improves achievable 
compression factors by an 

estimated one to three orders of magnitude and is well-suited to picture build-up or other iterative refinement. Combining these wavelet 
hlters and their inverses with previous theoretical work furthermore provides novel methods for calculating Fourier and other transforms 
In a preferred embodiment used to calculate the Fourier transform (1) and its inverse (34) applied to digital input data, the method replaces 
me fast Founer transform and its inverse and provides improvements in achievable accuracy. The new sampling method is inherently 
mu^tiscale and the invention thereby obviates the usual Nyquist constraim on the meaningful bandwidth in temis of the number of 
samples. Finally, the invention provides novel and efficient analog-to-ndigital and digital-to-analog interface. 
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METHODS OF DIGITAL FILTERING AND MULTI-DIMENSIONAL DATA COMPRESSION USING THE FAREY 
QUADRATURE AND ARITHMETIC, FAN, AND MODULAR WAVELETS 



It is well-known that digital filtering of digitally represented samples of analog data can 
be simplified and improved by performing the required processing in -some transformed 
domain. Using both a transform and its inverse to simplify the digital or other manip- 
ulation of data, one first transforms the data from its initial form to an intermediate 
form, then applies a specified manipulation, and finally applies the inverse transform to 
recover the data in its initial form from its intermediate form. By way of example and 
not by way of limitation: the initial form of the data may be digital samples of ana- 
log data, and the intermediate form of the data may be the coefficients of the Fourier 
transform; the initial form of the data may be the coefficients of the Fourier transform, 
and the intermediate form of the data may be the coefficients in some wavelet or other 
expansion; the initial form of the data may be an analog signal, and the intermediate 
form of the data may be some digital representation of it; the initial form of the data 
may be digital samples of an analog signal, and the intermediate form of the data may 
be a quantization of the digital samples. For instance, transform coding methods for 
data compression in signal or data processing are realized in practice as transform of 
data/quantization, storage or transmission, de-quantization/reconstruction. Another 
class of examples consists of finite-impulse response digital filters where the input 
data is filtered using the indirect calculation of convolution given by Fourier trans- 
form/multiplication/inverse Fourier transform. For analog data, sampling and filtering 
together determine fidelity and speed of manipulation. 

The invention disclosed here provides a new method for sampling analog or digital 
input data which may be used to calculate new wavelet transforms as well as their 
inverse reconstruction algorithms. This immediately provides novel and efficient meth- 
ods for data compression based on this new method of non-linear transform coding. 
In combination with these wavelet filters, the invention also provides new methods for 
calculating various classical transforms including the Fourier transform and its inverse. 
This immediately provides novel and efficient methods for digital filtering. By way 
of example and not of limitation, practical applications of the methods include non- 
speech audio compression, speech compression, speech recognition, speech synthesis, 
voice printing, audio filtering for hearing aids, still two-dimensional image compres- 
sion, moving video compression, video compression for purposes of telephony, precision 
Fourier analysis, precision trigonometry, denoising, interpolation, medical imaging, ge- 
ological imaging for recovery of oil or other resources, and other emission/detection 
apparatus. In specific applications, it may be necessary to store streams of input data 
in a buffer for separate processing or manipulation while concurrently collecting further 
input data whether analog or digital. It may also be useful in practice in particular 
applications to calculate and archive or otherwise store specific coefficients or constants 
and recover them from memory at run time rather than computing them at run time. 



In the landmark paper "Orthonormal bases of compactly supported wavelets, Com- 
munications in Pure and Applied Mathematics 41 , pp. 909-996 (1988), Daubechies 
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presented a class of \\ jiets which have been effectively applie- i numerous con- 
texts; see, for instance. US Patents 5,453,945 by Tucker et al., 5,6Ub,575 by Williams, 
as well as Daubechies' monograph "Ten Lectures on Wavelets" SIAM (1992) and its 
references. In another landmark paper, "An algorithm for the machine calculation of 
complex Fourier series", Mathematics of Computation 19, pp. 297-301 (1965), J. W. 
Cooley and J. W. Tukey described a simplified algorithm for the calculation of Fourier 
series on the computer. The utility of this fast Fourier transform FFT is well known 
with thousands of important applications and implementations; see, for instance, US 
Patents 3,871,577 by Avellar et aL, 3,881,100 by Works et aL, 4,051,357 by Bonnerot, 
and "The DFT, an owner's manual for the discrete Fourier transform", SIAM (1995) 
by W. Briggs and Van Emden Henson and its references. A requirement for both the 
FFT method and for the multiscale filtering of the Daubechies' wavelets is that the 
sampling of input data must be equally spaced though further methods of adaptive 
refinement of grids have also been developed in each context. 

Other examples of transforms whose efficient digital implementations are important 
and to which the invention is relevant include the transforms of Hilbert, Haar, Laplace, 
Bessel, Laguerre, Hermite, Chebyshev, Hotelling, Mersenne, and Fermat; see, for in- 
stance, US Patents 3,891,443 by Lynch et al. and 4,093,994 by Nussbaumer. 



Three new families of wavelets, respectively called arithmetic wavelets^ fan wavelets, and 
modular wavelets^ are defined on the standard circle of radius one in the complex plane. 
Any complex- valued function defined on this circle may be written essentially uniquely 
as a complex linear combination of each family, and the approximation of specified input 
data as a finite linear combination for each family leads to a corresponding wavelet filter. 
Each filter depends upon a new quadrature on the circle, called the Farey quadrature^ 
which relies on a novel method of non-equally-spaced sampling. The output of each 
wavelet filter plays the role in the invention of an intermediate representation between 
input data and any one of a number of useful output transformations of data which 
may be computed from the intermediate quantity. 

With either analog or digital input data, it is convenient to think of the given in- 
put data as spatial data. The inverse wavelet transform or reconstruction algorithm 
requires calculating spatial data from wavelet data and admits an especially conve- 
nient implementation: the values of the spatial function at certain non-equally-spaced 
grids of arbitrarily fine spacing may be computed exactly using a specified finite set 
of wavelet coefficients. This reconstruction algorithm combines with the wavelet fil- 
ter into a binary cascade giving an efficient procedure for data compression which is 
especially well-suited to discontinuous input data and to iterative refinement of out- 
put data. Furthermore, the methods extend directly to procedures for compression 
of multi-dimensional data. Computer coding for the reconstruction algorithm itself 
is sufficiently abbreviated and low-level that it might be transmitted together with 
compressed data, for instance, for down-loading from satellite to home computer. 

For the calculation of classical transforms of spatial input data, it is convenient to think 
of the desired output data, for instance, the Fourier coefficients, as frequency data. The 
transform from spatial to frequency data is accomplished in two main steps, the first 
of which approximates the input data by wavelets using the wavelet filter. The second 
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step is the calculation Jourier coefficients or other frequency dat . n:erra§"9f ^a^elets 
and requires formulas derived in previous theoretical work. Similarly, the calculation of 
the inverse transform depends upon the reconstruction algorithm and known formulas 
which express frequency data in terms of the wavelets. The calculation of a transform 
or its inverse admits an elegant implementation as a binary cascade. 



Brief Description of Figures 

FIG la illustrates the standard circle C of radius one in the complex plane, where 
C bounds the standard disk D of radius one. The complex numbers 4-1 and —1 
are also indicated, and the chord of C with these endpoints is labeled by the matrix 

FIG lb illustrates a chord in the top semi-circle of C labeled by a matrix A~ , 

together with two further ''descendant" chords of C labeled by the matrix products 

JJA and TA as indicated, where = and T = J^. The three chords 

determine a triangle inscribed in C, and the vertices of this triangle are indicated as 
complex numbers. 

FIG Ic illustrates a chord in the bottom semi-circle of C labeled by a matrix A = 
(c d)' together with two further "descendant" chords of C labeled by the ma- 
trix products V^A and T'^A as indicated, where = {^^^ l) "^"^ ~ 

1^)' three chords determine a triangle inscribed in C, and the vertices of 

this triangle are indicated as complex numbers. 

FIG 2 depicts a classical figure in mathematics called "the Farey tesselation in Klein's 
model of hyperbolic geometry" . It is constructed beginning from the chord in Figure la 
by recursively taking descendants as in Figure lb and Figure Ic, respectively, on the 
top and bottom semi-circle of C. The figure itself is classical, but the sampling of 
input data on C at the endpoints of the chords of the Farey tesselation in increasing 
generation is a novel aspect of the invention disclosed here. 

FIG 3a illustrates the circle C, the chord labeled / and the two triangles in the Farey 
tesselation on either side of this chord. The vertices of these triangles are indicated 
inside the circle as complex numbers. Outside the circle are given four explicit com- 
binations of cosine and sine functions, one such function next to each circular arc 
determined by the vertices, where 6 is the usual angular coordinate on the circle. A 
function t9/(5) is uniquely determined by the figure, where on each such circular arc, 
?9/(5) takes the values of the nearby combination of cosine and sine functions. The 
other parts Figures 3b-3e likewise determine analogous functions 'i9^(^) for other ma- 
trices ^ = ^d) '^^^^^ functions taken together form the system of wavelets on 
which the invention is based, and they are unique to this invention. 
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FIG 3b illustrates th- lord in the top semi-circle of C labeled th the matrix A, 
where a > c. As before, the figure determines a function on the circle denoted t9,4(0). 

FIG 3c illustrates the chord in the top semi-circle of C labeled with the matrix A, 
where c >= a. As before, the figure determines a function on the circle denoted 

FIG 3d illustrates the chord in the bottom semi-circle of C labeled with the matrix A, 
where a > — c. As before, the figure determines a function on the circle denoted t9>i(^). 

FIG 3e illustrates the chord in the bottom semi-circle of C labeled with the matrix 
A =, where — c >= a. As before, the figure determines a function on the circle denoted 

FIG 4 illustrates the notation near the chord labeled by the matrix A = 

the top semi-circle of C; this notation will be useful in several subsequent discussions. 

FIG 5 provides a flow chart for the main driving routine in calculating the Fourier 
transform. 

FIG 6 provides a flow chart for the main recursive routine in calculating the Fourier 
transform. 

FIG 7 provides a flow chart for the procedure of updating coeflacients in calculating 
the Fourier transform. 

FIG 8 provides a flow chart for the procedure of processing terminal arrows of the 
recursion in calculating the Fourier transform. 

FIG 9 provides a flow chart for the main driving routine in calculating the inverse 
Fourier transform. 

FIG 10 provides a flow chart for the main recursive routine in calculating the inverse 
Fourier transform. 

FIG 11 illustrates the region U consisting of all complex numbers with positive imag- 
inary part. Map the disk D toU via the function z i{z -h — 1); if two points in 
C are the endpoints of a chord in Klein's model of the Farey tesselation as in Figure 2, 
then draw the semi-circle in U passing through the corresponding two points which is 
perpendicular to the real axis. This produces another classical figure in mathematics, 
called the "Farey tesselation in the upper half-space model of hyperbolic geometry", 
which is indicated in Figure 11, 



Description of Preferred Embodiments 

In order to disclose the methods, it is necessary to develop some notation and termi- 
nology. To this end, consider the standard disk D of radius one about the origin in 
the complex plane with its usual complex coordinate 2 = x -t- iy, where i = \/— 1 and 
x,y are real numbers, together with its boundary circle C as illustrated in Figure la. 
An arrow is by definition an oriented chord of the circle C labeled by a two-by-two 

matrix ( ^ , ) , where the initial point of the chord is given by the complex number 
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€ C, and the final vint of the chord is given by eC, A licial arrow called 
the doe (short for "distinguished oriented edge'') is given by the diameter of C with 

endpoints -1,+1 labeled by the matrix ^ = ^ Figure la. Also required 

subsequently are the matrices 5 = ^ ^ U"^ ^ ^ ^ , = ^ ^ ^ ^ , for 

any integer n, which satisfy the identities = SUS = -T~^ STS = -C/"^ of 
matrix products. 

Inductively define two separate binary cascades, one of top arithmetic arrows and an- 
other of bottom arithmetic arrows^ as follows: 

Basis: The doe, which is illustrated in Figure la, is both a top arithmetic 
arrow and a bottom arithmetic arrow. 

Top Induction: As in Figure lb, given the top arithmetic arrow labeled by 
^ ~ d)' there are also top arithmetic arrows labeled by 

UA^i ^ ,^,) and r^=f^-^^ ^ + ^V 

Bottom Induction: As in Figure Ic, given the bottom arithmetic arrow 
labeled by ^ = ^ ^ ^ ^ , there are also bottom arithmetic arrows labeled 

It is a classical theorem in mathematics that the chords underlying any two arithmetic 
arrows meet only at their endpoints, if at all, and that the family of all such chords 
decomposes the disk D into an infinite collection of triangles with vertices in the circle 
C as is illustrated in Figure 2; this figure is called "the Farey tesseiation in Klein's 
model of hyperbolic geometry" . There is no suitable reference for this procedure in the 
literature, and it is for completeness that the expUcit recursive definition of arithmetic 
arrows has been given here. An arithmetic arrow is said to be of generation g if it arises 
firom g applications of the inductive clause starting from the doe in case of either top 
or bottom arrows; the doe has generation zero by convention. 

Consider the infinite set Q C C consisting of the endpoints of chords underlying all 
arithmetic arrows, and say a point of Q is of generation g if it is the endpoint of the 
chord underlying an arrow of generation g and no less; —1 and +1 are of generation 
zero by convention. There is a canonical enumeration of the generation g points of Q 
clockwise in C starting from —1 and the associated enumeration of Q itself by increasing 
generation. For example, the ordering begins with: 

-1 -4-1 -L- "1 - 2i 2-i 1 - 2i -1 - 3i 

, H-i, +2, 2, _j^22' -2 + 2' 2 + r i-f22' -i + sr 

Given input data on the circle C, the sampling of it at the points Q C C in this order 
of increasing generation will be referred to as the Farey quadrature. In practice, one 
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might interpolate given as required at the sampling points of tl "^arey quadrature 
or restrict the Farey quadrature to a circular arc in C. (In the subsequent discussion 
of mathematical basis, the sample points of the Farey quadrature will be seen to corre- 
spond to the rational points in the real line enumerated in increasing order, where the 
generation is the length of a corresponding continued fraction expansion.) 

Notice that the doe separates the triangle with vertices —1,-1-1,— i from the triangle 
with vertices —1, +1, -hi. In the same way, an arbitrary arithmetic arrow separates two 
triangles complementary to all the arithmetic arrows in D. The endpoints of chords 
of these triangles in various cases are illustrated in Figures 3a-3e as functions of the 

matrix ^ = labeling the arithmetic arrow. Figure 3a depicts the doe itself 

with j4 = /, Figure 3b and 3c depict top arrows with a > c and c > a, respectively, and 
Figure 3d and 3e depict bottom arrows with a > — c and — c > a, respectively. The 
notation inside the circles indicates the endpoints of the chords as complex numbers, 
and the notation outside the circles in each case of Figures 3a-3e is y^t to be explained. 

To this end, define a trigonometric function to be a real-valued function 

f{0) = a cos e -h /3 sin e H- 7, 

where 9 is the usual angular coordinate on the unit circle C and a, 7 are real numbers. 
Next, define a piecewise trigonometric function on the circle to be a function f(6) so 
that there is a finite collection of circular arcs /j C C, for = 1, . . . , J, any two of which 
meet at common endpoints in C, if at all, together with a collection of trigonometric 
functions tj(9) for j = 1, . . . , J, so that f{e) = tj(9) if 9 e Ij. In other words, C is 
decomposed into a finite number of non-overlapping circular arcs, and on each such 
arc, / takes the values of a trigonometric function. 

There is a particular family of piecewise trigonometric functions, one such arithmetic 
wavelet i9yi(5) for each matrix A labeling an arithmetic arrow. These functions are 
defined in Figures 3a-3e, where the endpoints of the circular arcs are as described 
before, and each trigonometric function is written outside the circle next to its corre- 
sponding circulsur arc. There are thus five cases in the definition of arithmetic wavelet 
corresponding to the five cases in Figures 3a-3e already discussed, and in each case, the 
arithmetic wavelet itself is a piecewise trigonometric function defined using exactly four 
circular arcs. Arithmetic wavelets are a fundamental new ingredient of this invention 
and no doubt seem entirely ad hoc and unmotivated at this juncture of the discussion. 
The mathematical basis for them will be given later, but it is worth pointing out now 
that each arithmetic wavelet '0Ai9) is a once-continuously diff'erentiable function taking 
value zero at the points -1,4-1, -i of the circle C, that is, for 0 = 0, -7r/2, -tt, except 
for '^^-l{9) and ^7»-i(5) which take the value zero only at the points ±1. 

In order to explicate the definition of arithmetic wavelets, notice that any two-by-two 
matrix A= determines a self-mapping Ma : C C ol the circle C as follows: 

{ta + 6) - i{tc-^d) 



Ma 



{ta-¥b) + i{tc-\- d)' 



and in fact, the image of the endpoints, -1 and -i-1, of the doe under Ma are the 
endpoints of the arithmetic arrow with label A. Construct from the function dj{0) 
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defined on C in Figure L .he corresponding vector field on C given ^liO) where 
^ is the usual unit tangent vector field to C 

There is an explicit formula to be given presently for the transformation of suitable 
vector fields under the change of coordinates on C given by Ma- Namely, the vector 
field t{e)'§Q = [a: cos 0 + ^ sin ^ + 7)^ transforms under Ma to the vector field 
where 

i [2{cd - ab)P + (d^ - b^)(a + 7) + (a^ - c'')(a - 7)] cos 0 
= I -h [(ad + 6c)/? -f M(a + 7) - ac(a - 7)] sin <9 |. 
+ ~ [2{cd + ab)p + + b^){a + 7) - (a^ + c2)(a - 7)] 

More generally, given a piecewise trigonometric function f(0) taking the values of the 
trigonometric function tj{6) on the circular arc Ij as before, the corresponding vector 
field f{9)^ transforms under Ma to the vector field which on the circular arc MA{Ij) 
is given by tf{9)^. 

Changing coordinates on C by Ma using this formula transforms the vector field 
^/(^) ^ to another vector field dAiO)^, thereby defining the function on C 

Finally except for A = i7~^ and A = T~'^ , the corresponding arithmetic wavelet is 
given by 

^a{0) = dA{0) - [or cos e + /? sin 6 + 7], 
where a,P,y are chosen so that 79(7r) = d{0) = d(—7r/2) = 0; the two special cases are 

du^i{9) = i9t/-i {6) -f- 2 sin e and (9) = ^r-^ W - 2 sin 6> 

(which are defined in this manner to guarantee a symmetry of the upper and lower 
semi-circles of C). This is the abstract definition of arithmetic wavelets which is made 
explicit in Figure 3. 

It is furthermore convenient to define 

(dAi9y, ifyl^^/-^T-^ 

['dT-^(9); ifA = T'\ 
so that for every label A on an arithmetic arrow, 7?>i(7r) = ??>i(0) = 79^(— 7r/2) = 0. 

The arithmetic wavelets enjoy an important renormalization property, namely, 

^ba{Ma{9)) - MM^b{9). 

where the prime denotes the usual derivative and where it is required that either both 
matrices B and A are matrix products of U'^^,T'^^ or both matrices B and A are 
matrix products of W^^T"^. 

Just as a function / ~ f{9) may be expressed in its representation as a Fourier series / ^ 
J2n ^Vifi'"^? SO too it may be expressed as a series f ^ eA'^A{9) of arithmetic wavelets, 
which together form a linearly independent set. The sum is over all arithmetic arrows, 
and the coefficients ca are called the arithmetic wavelet coefficients. The calculation 
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of arithmetic wavelet c ncients from the input data / is called th rithmetic wavelet 
filter. The arithmetic wavelet coefficients are not quite uniquel3'' determined by / as 
will be explained. 

The reconstruction algorithm or inverse wavelet transform provides a method for cal- 
culating function values at points of the circle from arithmetic wavelet coefficients. To 
determine the value of the sum ^a'^a at a specified point g G Q Q draw a line 
segment from g to the origin in D; this segment is contained in a finite collection of tri- 
angles complementary to all arithmetic arrows, whose boundary edges taken together 
form a finite set of chords; the corresponding finite set of arithmetic arrows A enu- 
merates the only coefficients eA which affect the value of YIa^s ^a'^a at g. This is an 
important finiteness property of arithmetic wavelets which is crucial to the efficiency of 
the reconstruction algorithm. 

This completes the general discussion of arithmetic wavelets, the first family of new 
wavelets disclosed here, and the two further families of new wavelets are next defined. 

For each label A on an arithmetic arrow, there are two fan wavelets 

4>A{e) = J2^ur^^{0h 

n>0 

<PsA(O) = J2^T-nA{0), 

as well as two modular wavelets 

n>0 



n>0 

It follows immediately from the definitions that there are the following relations among 
the wavelets for any label A on an arithmetic arrow. 

^a(0) - ^u-^A{0)-Md). .4>sa{0) = i^u-^sA{0)-i^sA{e), 
^A{e) = ct>A{e)-<f>uA{e) = <f>sA{0) - 4>asA{e) 

= ^U~^A{0)-2lPA{0)^1pUAi0) = ^U-^SAi0)~'^i^SAi0)-^i^USAi0)- 

On the other hand as is not at all obvious from the definitions, cpj and 'tpi are given by 
the following explicit formulas as piecewise trigonometric functions 

{2 sin 0; if 0 < 0 < tt, 

-2 cos 0 - 2 sin 5 - 2; if n < 0 < 37r/2, 
-2 cos e -I- 2 sin 5 + 2; if 37r/2 < 0 < 27r, 

and 

, f 2-2 cos 0; if 0 < 5 < TT, 

^^(^)-=^(o; ifnl0l2.. 

In fact, modular wavelets arise from ipj in exactly the same manner that arithmetic 
wavelets arise from t?/, namely, ^/(5)^ transforms under Ma to the vector field 
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"PAWae^ ^a{0) — ,(^) - [a cos 5 + /? sin ^ + 7], where a, 7 are chosen to 
guarantee that ^^^(0) = i^'a(^) = 7i-u(~7r/2) = 0: transforming using the 

change of coordinates given by Ms a likewise gives rise to iPsa(0) - The fan wavelets 
also arise from (pj in the analogous manner. 

In contrast to arithmetic wavelets, neither the fan nor modular wavelets are linearly 
independent, and in each case, there is one linear relation for each arithmetic arrow, 
namely, 

4>A — <i>UA = <f>SA — <t>USA^ 
'IpU-^A - ^i^A + TpUA = IpU-^SA - ^'fpSA + "ipUSA- 

On the other hand, the following collection of modular wavelets forms a basis 
{ipA ' A labels a top arithmetic arrow and c> a} 
U {'ipsA • A labels a top arithemtic arrow and a > c} 
U {ipA ' A labels a bottom arithemtic arrow and — c > a} 
^ {^5>i : A labels a bottom arithemtic arrow and a > — c}; 

there is an analogous basis for the fan wavelets. 

Both fan and modular wavelets enjoy precisely the same finiteness property as arith- 
metic wavelets. Furthermore, both fan and modular wavelets enjoy renormaUzation 
properties, namely, 

<t>BA{MA{e)) = M'^{e) (f>Bi0). 

where the requirements on A,B are as before. 

This completes the necessary abstract definitions. The methods will first be presented 
using arithmetic wavelets for real-valued input functions with extensions to modular 
and fan wavelets and complex- valued or multi-dimensional input/output data to be 
described later. As a point of notation throughout these discussions, there will be the 

tacit assumption that the entries in the matrix denoted A are given by 




The Fourier Transform 

Included in the section entitled "Source Codes" is an implementation in the computer 
language C of a preferred embodiment of the method described in this section employing 
arithmetic wavelets. 

Suppose that f(6) is some specified real-valued input function defined on the circle. 
The calculation of the Fourier transform of / proceeds by two main steps: 

Stepl : Choose a finite set S of arithmetic arrows and approximate / ^ 
5Da€5 ^a'^a calculating the coefficients ca in a manner to be described in 
terms of the values of / at specified points using the Farey quadrature. 

Step 2: Substitute into the expression in Step 1 the known Fourier expan- 
sions ^9/1 (5) ~ Yin s*"^ in order to derive the approximation 

Cn = 5^ c^f. 
Aes 



9 



wo 00/441104 

^ uu/*(^^iiu^ PCTAJS99/30584 



The coefficients require a Step 2 are known theoretically and gi\ for n^^ 0,±1 by 

1(0 — a) + i{a — c) 
+ 2(c2 + d2)[^1" +.2(a2 + fc2)r^l" 

-|(. + a)»■^(6 + .flf^i±§^4^i±41" 

+ +z(a + c)J 

Furthermore, the coefRcients CQ,cf,c^^ are given bj' 

TTC^i = 0h(c2 -S± 2icd) + er[(6 - d)2 - (a - c)^ if 2i(6 - d){a - c))/2 
+ Otio? -b"^ ±2iab) + 6t[{b + d)'^ -{a + cf^:2i{b + d){a + c)\/2, 

wc^ = 29h(c^ + d'') - 0,[(b - df + {a - cf] 
+ 20tia^ + b'^) - 9i[{b + d)^ + {a + c)% 

where the angles are 

It may be useful in practice in particular applications to calculate and archive or oth- 
erwise store the constants required in the previous expressions for and recover them 
from memory at run time rather than computing them at run time. 

The Fourier coefficients co,c±i of / thus require special treatment which does not 
present additional challenges. For simphcity, suppose that the input function f(0) 
has these three coefficients fixed by the condition that f{0) = 0 for 0 = 0, -7r/2, -tt. 
More precisely, if f{0) is any real- valued function on the circle, then let J (9) denote its 
normalization defined by f {9) = f{e) + a cos ^ + ^ sin 5 + 7, where a, jp, 7 are chosen 
so that / takes value zero at 0, -7r/2, -tt. The specified simplifying assumption on / 
thus amounts to the condition that / = /. 

The procedure for calculating wavelet coefficients is shghtly different depending upon 
whether the corresponding arithmetic arrow is a top or bottom arrow. However, the 
antipodal map 0 ^ it + 0 transforms the bottom of the circle C to the top in such a 
way that it suffices to discuss only the case of top arrows. More precisely, the algorithm 
for bottom arrows arises from that for top arrows by: 

o replacing the input data value f{6) with the input data value f{9 + tt), 

o replacing the matrix A in the described manipulations for top arrows 
by the conjugate matrix SAS for bottom arrows, and 

o replacing the trigonometric function Qfcos 0 -h /?sin 0 + 7 fo^^ top arrows 
by the trigonometric function — acos 9 - /3sin 5 -f 7 for bottom arrows. 
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riie algorithm and cor. iter coding for bottom arrows is therefi. derived without 
difficulty- from that for top arrows, and the subsequent discussion is specialized without 
loss of generality to top arrows. 

In addition to the function f{d) defined on the circle, suppose that there is also specified 
a bandwidth of interest N (i.e. the method should calculate the Fourier coefficients Cn 
only in the range \n\ < N) and a tolerance EPS (where contributions to the specified 
Fourier coefficients are regarded as negligible if they are of magnitude less than EPS). 
There are numerous variants of the method which depend upon other parameters which 
may be tailored to a given application, and some these variants will be discussed further. 

Steps 1 and 2 are merged into a single binary cascade which keeps track of an ongoing 
approximation to the Fourier coefficients as follows. A cascade element is called an 
arrow-structure and is defined to consist of the specification of 

o an arithmetic arrow A of generation 

o the wavelet coefficient ca on '3 a in the approximation of /, 

o the coefficients a, /?, 7 of a trigonometric function defined as follows. The end- 
points of the chord corresponding to A decompose C into two circular arcs, ex- 
actly one of which contains none of the points —1, -fl, -2. The sum JZbj^a ^b'&b 
over all arithmetic arrows B of generation at most g except for A itself is given 
on this interval by the particular trigonometric function a cos ^4-/3 sin ^ -f 7. 

The coefficients a,/3,7 in an arrow-structure keep a lagged/updated running tally of 
the overall effect of what has come before it in the cascade; as a result of this technique, 
the method requires essentially no memory other than the storage of the running ap- 
proximation to Fourier coefficients and the stack required for the recursive computation 
in the cascade of arrow-structures. 

The next technical point involves a regularization in the calculation of arithmetic 
wavelet coefficients; this regularization is required because the approximation to / as a 
linear combination YIa^s ^a^a is not uniquely determined. Given an arrow-structure 
in the specified notation and referring to Figure 4, when sampling the input data at 
the point in the Farey enumeration with corresponding angular coordi- 

nate 5o, the one real input datum f{0o) must determine the two coefficients euA and 
e-TA^ According to the formulas presented in Figures 3a-3e, these new coefficients are 
constrained by 



KOo) = a 



jb + d)"^ - (q + c)^ 
_(6-hd)2-f-(a-hc)2 



2(a + c) (6 + d) 



(6-f d)2+(a-hc)2_ 
MeuA - cta) 8 ca sgn (a - c) 



H-7 



where sgn (a - c) is a sign defined to be +1 if a > c and -1 otherwise. There remains, 
however, one real degree of freedom in the specification of euA, ^ta, and it is eliminated 
by demanding the best least-squares fit to the values of / at the next generation of points 

(26+d)4-i(2a+c) (64-2d)+t(a-t-2 c) ^. , j ^ ^ n .„ , 

(26+rf)-i(2a-(-c)' (6+2d)-t(a4-2c) respective angular coordmates 61, $2, as illustrated 

in Figure 4. Updating the lagged trigonometric coefficients in accordance with the 
formulas presented in Figures 3a-3e produces trigonometric functions Vi(^) and ¥2(6) 
which give the values of the ongoing approximations at 0i and 82, respectively, and 
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whose coefficients are e icit inhomogeneous linear functions of e; bta, namely, 

ViiOi) = vie + V2U euA + vit er^, for i = 1,2. 

The coefficients are explicitly described most concisely in pseudo-code, where e = e^, 
alp = a, bet = /5, gam = 7, bpd = b + d and ape = a + c: 

if (a > c){ 

axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 
ayl=2*(b*b-a*a); 
bxl=bet-f-(c*d-a*d-b*c-a*b)*4*e; 
byl=4*a*b; 

cxl=gam+(2*(a*c-)-b*d)-c*c-d*d+a*a+b*b)*2*e; 

cyl=-2*(a*a+b*b); 

ax2=alp-f-4*(d*d-c*c)*e; 

ay2=2*(c*c-d*d); 

bx2=bet+8*e*c*d; 

by2=-4*c*d; 

cx2=gam-4*e*(c*c+d*d); 

cy2=2*(c*c+d*d); 

} 

else{ 

axl=alp+4*e*(a*a-b*b); 
ayl=2*(b*b-a*a); 
bxl=bet-8*e*a*b; 
byl=4*a*b; 

cxl=gam+4*e*(a*a4-b*b); 
cyl=-2*(a*a-fb*b); 

ax2=alp-|-(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c-l-c*d)*4*e; 
by2=-4*c*d; 

cx2=gam+(a*a-|-b*b-2*(a*c+b*d)-c*c-d*d)*2*e; 

cy2=2*(c*c+d*d); 

} 

vlc=cxl+ (((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl 

-l-2*(b+bpd)*(a-l-apc)*bxl)/((b+bpd)*(b+bpd)-l-(a-fapc)*(a+apc)); 

vlt=cyl-{-(((b+bpd)*(b+bpd)-(a-l-apc)*(a-l-apc))*ayl 

-t-2*byl*(b-l-bpd)*(a-l-apc))/((b+bpd)*(b-l-bpd)4-(a-l-apc)*(a+apc)); 

vlu=8./((b+bpd)*(b+bpd)-|-(a-l-apc)*(a4-apc)); 

v2c=cx2-|-(((d+bpd)*(d+bpd)-(c+apc)*(c-l-apc))*ax2 

+2*(d+bpd)*(c+apc)*bx2)/((d-|-bpd)*(d+bpd)-l-(c+apc)*(c+apc)); 

v2u=cy24-(((d4-bpd)*(d+bpd)-(c+apc)*(c-l-apc))*ay2 

+2*(d-|-bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c-)-apc)*(c+apc)); 

v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c-l-apc)); 

It is easy to minimize [f{9i) - Vi{ei)f + [f{02) - ¥2(02)]^ as a function of euA,eTA 
subject to the given constraint, and this procedure provides a suitable regularization. 
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Weighted least-squares .s or other statistical treatments using i. it data values at 
further generations allow aiternati\-e regularizations. 

Two other simple regularizations which have also been useful are: to require that 
euA = 0 if c > a and cta = 0 if a > c, or, to require that euA/^TA = —[a{a + c)-\-b{b + 
d)]/[c(a-fc) + d(6 + c?)]; the latter reguiarization is related to imposing differentiability 
of the approximation at Oq- 

The final technical point involves the stopping criteria and terminal processing for the 
cascade of arrow-structures. The basic stopping parameter is NVAN, where a branch 
of the cascade terminates whenever there have been NVAN consecutive generations 
of offspring whose contributions to all coefficients Cn in the bandwidth N have been 
of magnitude at most EPS. There are further parameters governing stopping criteria 
which may be introduced depending upon the scale and resolution of features in the 
input data, including a minimum or maximum number of generations, the size of the 
angle subtended by the chord of a terminal arrow-structure, and so on. 

Another useful stopping criterion is to require that |e^| W'OaW be negligible for sev- 
eral generations, where W^aW denotes the L^-norm of 3 a for which there is the a 
priori estimate W^aW < 10|ac -h that is, one terminates the cascade when 

NVAN consecutive generations of coefficients ca have satisfied the condition that 



When NVAN generations of offspring contribute negligibly to the specified bandwidth, 
the cascade is terminated with those offspring of greatest generation. Suppose that A 
labels the arrow of such a terminal offspring, adopt the notation in Figure 4, and let 
J C C denote the circular arc with endpoints ^37^, ^rff which contains none of the 
points —1, —2, -f 1 € C. The final update n = alpz cos 6 + beti sin 0 4- gami, for i = 
1,2, of the running trigonometric function on J in accordance with Figures 3a-3e is 
prescribed in pseudo-code with notation as before as follows: 



if(a > c){ 

alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); 

bet l=bet+e*4* (c* (d-b)-a* (d+b)) ; 

gam 1 =gam+e*2* (a* (a+c) +c* (a-c) +b* (b-l-d)-d* (d-b) ) ; 

alp2=alp+e*4* (d-c) *(d-f-c) ; 

bet2=bet+e*8*c*d; 

gajn2=gam-e*4*(c*c-l-d*d); 

} 

else{ 

alpl=alp-l-e*4*(a-b)*(a+b); 

betl=bet-e*8*a*b; 

gam 1 =gam-l-e*4 * ( a* a-l- b*b) ; 

alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b-l-d)); 
bet2=bet+e*4* (a* (d-b) -1-c* (b-l-d) ) ; 
gam2=gam-He*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b-|-d)); 
} 

Now serially letting t denote the trigonometric functions ri and T2, again adopt the 
notation of Figure 4, where A labels the arrow corresponding to r. The updated r 
is compared with another trigonometric function a defined by extrapolation which is 



^ EPS |ac + 6d|5/4. 
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uniquely determined b;, .^manding that it agree with / at the thr^ .ngles comprising 
the next two generations, B = ^0,^1,^2 in the notation of Figure 4. The Fourier 
coefficients of the piecewise trigonometric function which takes values a{9) — r{0) on 
the circular arc J and value zero otherwise are easily computed using standard formulas 
and added to those of the ongoing approximation. 

Other algebraic or statistical samplings of higher generation offspring allow alternative 
final processings. 

For instance, one might save computational expense and simply terminate the cascade 
with no final processing at all. 

Furthermore, the terminal arrow-structures could be stored for restart or iterative 
refinement capabilities. 

One might alternatively record in an arrow-structure different transforms of the trigono- 
metric function defined by a,/3, 7. For instance, it is useful to take advantage of the 
renormalization property of wavelets and change these coefficients by transforming the 
vector field [a cos ^-1-/3 sin 64-7]^ by the change of coordinates M^^ This technique 
avoids the necessity of calculating certain large integers attendant to the calculation of 
wavelet coefficients as will be illustrated later. 

One favorable aspect of the method is its advantageous mix of floating-point and integer 
operation types. Moreover, except for the coefficients Co, c±i, the implementation of the 
algorithm is purely algebraic, that is, requires only addition and multiplication. Fur- 
thermore, by its very nature as a cascade, the algorithm is amenable to parallelization 
and efficient hardware implementation. 

This completes the description of the implementation of the method of calculating 
Fourier coefficients using arithmetic wavelets. 

Several points will next be addressed which are special to the implementation of the 
methods using modular wavelets. Again there is a basic Step 1, the calculation of 
the modular wavelet transform where the input function is approximated as f{6) ^ 
X^/4pi4^>i(^), and a beisic Step 2, substitution of known expressions for the Fourier 
coefficients of the modular wavelets; these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. 

To elucidate Step 1 for modular wavelets, refer again to Figure 4 for the notation near 
the arithmetic arrow labeled by the matrix A, The single value of the input func- 
tion /(^o) together with the ongoing calculation of the updated trigonometric function 
Or cos 5 4- /? sin & -i- 7 in the arrow-structvue this time uniquely determines the modular 
wavelet coefficient 

^ f{0o) - {oc cos 90-^13 sin Bp + 7) 

since tpsiOo) = 0 for every labeling matrix B of generation greater than A] there is 
thus no need for regularization with modular wavelets. 

As to Step 2, recall that a basis of modular wavelets was given before by a collection 
^b(^) of wavelets where A is the label on an arithmetic arrow and either B =: ^4 or 
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B = SA, In each case, 3 can calculate that tpsid) has Fourier e: nsion 



where 



|n|>l 




<' b-m \" _ / tf-tc N" 
\b+iaJ \d+icJ 



; if 5 = ^1, 

; if B = 



,B ^ / -4i ^ (Irif ifB = A 



and 

° l-4f^-2(a2+62)[tan-^j2s|^-tan-i5l£i,]; if B = 
2nr = 2;rd^ - / id)2[tan-^^ - tan-^^lf^]; if B = A, 



Thus, the method of calculating Fourier transforms using modular wavelets is easily 
derived from the method using arithmetic wavelets. In fact, it is a simple matter to in- 
corporate the two minor modifications just described and alter the included source code 
employing arithmetic wavelets to produce source code employing modular wavelets. 

In the same way for fan wavelets, there is again a Step 1 and Step 2, which are merged 
into a binary cascade of arrow-structures. As for arithmetic wavelets, Step 1 for fan 
wavelets again requires a regularization scheme. As to Step 2, the formulas for Fourier 
coefficients of fan wavelets can be derived from those given before for modular wavelets 
using the identity 0^(5) = ipAiO) - i^UA{6) mentioned before. Again, the included 
source code is easily modified to employ fan wavelets. 



Flow Charts for the Fourier Transform 

Flow charts for various procedures used in the calculation of Fourier transforms are 
presented in Figures 5-8. For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the renormalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are referred to in these figures simply as "arrows" . 

In Figure 5 is presented a flow chart for the main driving routine. The input data is 
normalized as indicated in program segment 1. There are separate recursions estab- 
lished in program segments 2 and 3 for the top and bottom of the circle respectively. 
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Each recursion is estab' .ed with a call to the subroutine gener, w' h then calls itself 
in turn. The Fourier coefficients are stored internally with an overall suppressed factor 
of TT, and output data normalization of multiplying by tt is accomplished in program 
segment 4. Output Fourier coefficients are finally displayed before exiting in program 
segment 5. Program segments 2 and 3 are entirely independent and could be performed 
in parallel; more generally, each of program segments 2 and 3 could be decomposed 
further into multiple parallel procedures. 

A flow chart for the main recursive routine gener is given in Figure 6. The procedure 
starts with a test in program segment 6 to determine if: 

o the argument envy >= NVAN, in which case there have been NVAN 
consecutive generations of negligible contributions, and 

o the argument generation >= MING, in which case there have been a 
required minimum number of iterations in the recursion. 

If both of these inequaUties hold, then the procedure passes to program segment 15, 
where the cascade is terminated and the output data corrected with a call to the 
subroutine prune. 

In the contrary case that the recursion should continue, the procedure passes to program 
segment 7, where the descendant arrows in the cascade are determined using the least- 
squares fit to the next generation of data as described before to compute the regularized 
wavelet coefficients of the descendants. These are combined with integer calculations 
to update the lagged trigonometric functions. The procedure then passes to program 
segment 8, where the ongoing approximations to Fourier coefficients are updated to 
include contributions from the two descendant arrows with a call of the subroutine 
Charles for each descendeint. The procedure continues with a test in program segment 
9. If the contribution calculated in the subroutine charles for either descendant to any 
Fourier coefficient in the bandwidth was non-negligible, then the recursive argument 
envy is set to zero in program segment 10, In the contrary case that both contributions 
to all Fourier coefficients in the bandwidth were negligible, the control parameter envy 
is increased by one in program segment 13. 

In any case, the procedure passes to program segment 11, where there is a test on the 
number of generations. In case too many iterations of the recursion have occurred, 
it is terminated in program segment 12, and suitable warning of non-convergence of 
the method is written to the output. In the contrary case that the maximum number 
of generations has not been reached, program segment 14 establishes the recursion by 
calling gener once for each of the descendant arrows with the updated control parameter 
envy and an incremented generation. 

In Figure 7 is presented a flow chart for the subroutine charles. Calling the subroutine 
Charles has the effect of updating the ongoing approximations to the Fourier coefficients 
for a single argument arrow and returning a flag which keeps track of whether any 
such contribution has been non-negligible. The flag is cleared in program segment 16, 
and several preliminary calculations are accomplished in program segment 17. The 
procedure passes to program segment 18, where it is determined whether to calculate 
the 0,-i-l,-l Fourier coefficients. If these are to be calculated, then the procedure 
passes to program segment 19, which calls a subroutine to update these three Fourier 
coefficients. Program segments 20 and 21 set the flag as required depending upon 
whether these three contributions are non-negligible. 
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In any case, the proce re passes to program segment 22, If the op over the band- 
width is complete, then all Fourier coefficients have been updated, so the subroutine 
Charles returns the flag in program segment 23. In the contrary case, the current Fourier 
coefficient is updated in program segment 24. The flag is set as required depending on 
whether the current contribution was non-negligible in program segments 25 and 26. 
The iteration over bandwidth is established by the update in program segment 27 and 
the return to program segment 22. 

In Figure 8 is presented a flow chart for the subroutine prune, which modifies the ar- 
ray of Fourier coefficients when terminating the cascade. The final update of lagged 
trigonometric functions to produce T2 is performed in program segment 28, There is 
one such function for each possible descendant of the argument arrow. The procedure 
passes to program segment 29, where further input data is sampled in order to com- 
pute two trigonometric extrapolations guo^, one such extrapolation for each possible 
descendant of the argument arrow. Each function - r,, for i = 1, 2, is truncated as 
described before, and the Fourier coefficients of the truncated functions are added to 
the ongoing approximations of Fourier coefficients in program segment 30: Program 
segments 31 and 32 implement the calculation of 0,-1-1,-1 Fourier coefficients if desired, 
and the subroutine prune is terminated with the return in program segment 33. 



Included in the section entitled "Source Codes" is an implementation in the computer 
language C of a preferred embodiment of the method described in this section employing 
arithmetic wavelets. 

Since the method for calculating the inverse Fourier transform is similar in spirit and 
execution to that for the Fourier transform described in detail before, it need only be 
briefly discussed. The input data includes the specification of a collection of Fourier 
coefficients in a given bandwidth iV. 

There are again two main steps in the calculation: 

Stepl : Calculate arithmetic wavelet coefficients from the given Fourier 
coefficients Cn- 

Step 2: Use the wavelet coefficients from Step 1 and the reconstruction al- 
gorithm to output values of the function YLa ^/I^a(^) at the Farey quadra- 
ture points on the circle. 

For Step 1, there is again an explicit formula from previous theoretical work. Deep 
mathematical results prove that for n 0, ±1, 

e'"^ = cos nO i sin nd 



where the sum is over all arithmetic arrows, the underlying chord of which has endpoints 
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^, 77 € C, and 



Co = 



-1, 


n 


= 0(4) 


0, 


n 


^1(4) 


-1, 


n 


= 2(4) 


0, 


n 


= 3(4) 



Ci = 



0, 


n 


= 0(4) 


-1, 


n 


= 1(4) 




n 


= 2(4) 


0, 


n 


= 3(4) 



c-1 = 




n = 0(4); 
n = l(4); 
n = 2(4); 
n ~ 3(4). 



Substituting this expression for e^"^ into the given Fourier expansion f{d) ^ J2 



n 



inO 



gives 



where 



and 



A 



Ce = Ci- ^ Cn cP, for ^ 



0,±1. 



This expression /(<9) « + cie^^ + clje"*'^ + 79^(5) of /(/9) constitutes Step 1 

in the calculation of inverse Fourier transform. 

In fact, the calculation of wavelet coefficients in Step 1 can be implemented using 
mostly integer operations to improve run-time dramatically. To see this, notice that 
the formulas ^ = 7] = for the endpoints of the chord underlying the arrow 

labeled by the matrix ^ = may be substituted into the expression for ca to 

yield 

where 

n ^ [(bd-ac) +i(arf + bc)]^ 

^n{{bd -hac-h iT + {bd + ac - t)"} -f i{bd -h ac) + ac + z)" - (fed + ac - i)^}J , 

In the obvious implementation of this formula to calculate e% using integer operations, 
the intermediary integer expressions grow too large to be practical for n large if the 
arrow corresponding to A is of large generation. 

This difficulty is overcome by wavelet renormalization, where these attendant large 
integer calculations are obviated entirely as is illustrated in the included source code. 

Step 2 in the method for calculating inverse Fourier transforms depends upon the 
reconstruction algorithm to output function values j>n the circle taken by the linear 
combination /(i9) ^ ^ c\e'^ + clic"*^ + YLa OAiO). The finiteness property of 
the reconstruction algorithm allows the exact calculation of these function values at 
the sample points Q C. C \n their ordering determined by the Farey quadrature. The 
two steps are again conveniently merged into a single binary cascade as follows. 
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There is only one con 1 parameter, which is called SCALE, v. :e the method is 
required to determine at least one output value in each circular arc subtending an 
angle SCALE. Thus, as SCALE is decreased the grid of output values on the circle is 
refined. A binary cascade is established using arrow-structures which is similar to that 
for the Fourier transform before. However, the wavelet coefficients are given in this case 
in closed form by the formula for ca (so no least-squares fit is required). Because of the 
finiteness property of the reconstruction algorithm, exact output values are determined 
at each stage of the recursion. The cascade is terminated with an arrow whenever the 
chords underlying both descendant arrows subtend angles less than SCALE. Output 
data can be written to an array either when calculated to store data in the ordering of 
the Farey quadrature or when terminating the cascade to store data in the standard 
ordering on the circle. Furthermore, the terminal arrow-structures could be stored for 
restart or iterative refinement capabilities. 

This completes the description of the implementation of the method of calculating the 
inverse Fourier transform using arithmetic wavelets. 

Several points will next be addressed which allow the extension of the methods to 
modular and fan wavelets. Again there is a basic Step 1, the calculation of wavelet 
coefficients from Fourier coefficients, and a Step 2, the reconstruction algorithm; in 
each case of fan or modular wavelets, these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. The finiteness condition on fan 
or modular wavelets mentioned before renders Step 2 entirely analogous to that for 
arithmetic wavelets. As to Step 1, simply substitute the identities '^a = (f>A- <t>UA = 
'ipUA - 2i}A + i^u-^A described before into the expression above for e*'*^ in terms of 
arithmetic wavelets to immediately derive corresponding expressions in terms of fan 
and modular wavelets. 

Thus, the method of calculating inverse Fourier transforms using fan or modular 
wavelets is easily derived from the method using arithmetic wavelets. In fact, it is 
a simple matter to incorporate the minor modifications just described and alter the 
included source code employing arithmetic wavelets to produce source code employing 
fan or modular wavelets. 



Flow charts for various procedures used in the calculation of inverse Fourier transforms 
are presented in Figures 9-10. For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the re-normalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are referred to in these figures simply as "arrows". 

In Figure 9 is given a flow chart for the main driving routine. Modified Fourier coef- 
ficients Co,Ci,ci.i are computed in program segment 34. The formula given in Step 1 
for the inverse Fourier transform is applied in program segment 35 to calculate the 
arithmetic wavelet coefficients in terms of Fourier coefficients for a particular family 
of nine arrows. These nine arithmetic wavelet coefficients are required to initiahze the 
trigonometric functions for recursions established in program segment 36 with calls 
to the recursive subroutine gener. Recursions with initial conditions corresponding to 
^ — U,T,T-^, T~^T"^Z7~\ are undertaken in program segment 36, These 
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recursions, as well as fi ler possible subrecursions derived from m, are indepen- 
dent and may be undertaken in parallel. This elaborate driving routine is convenient 
because e*"^ is expressed in Step 1 ofthe inverse Fourier transform method in terms 
of the normalized arithmetic wavelets t^a which differ from the arithmetic wavelets i?^ 
themselves only for A = U~^^T'~^, 

The procedure passes to program segment 37 to perform an overall correction of the 
function values using the modified Fourier coefficients Co,ci,c'_i. Output data is dis- 
played and the procedure terminated in program segment 38. 

In Figure 10 is given a flow chart for the main recursive routine gener. In the cas- 
cade, the two descendant arrows of the argument arrow are generated in program 
segment 39 employing Step 1 of the inverse Fourier transform method to calculate the 
wavelet coefficients of the descendants. These expressions are then used to update the 
lagged trigonometric functions in program segment 39. A test is performed in program 
segment 40 to determine if the chord underlying each descendant arrow subtends a 
sufficiently small angle, as estimated by the comparison of an integer quantity with 
SCAT; the relation with the control parameter SCALE discussed before is given by 
SCAT = [expsinh-^(l/SCALE)]2, In case one of the chords subtends too large an 
angle, the procedure passes to program segment 41. The recursion is established in 
program segment 41 with two calls to gener, where the arguments are given by the 
two current descendant arrows. In the contrary case, the procedure performs the final 
update of the lagged trigonometric functions in program segment 42, then stuflFs the 
output array with the new function values in program segment 43, and finally returns 
in program segment 44. 



Data Compression 

As with any method for data compression based on transform coding, five basic ma- 
nipulation are required, namely: 

o wavelet filtering, 
o quantization, 

o storage and retrieval or transmission 
o de-quantization, and 
o reconstruction. 

The wavelet filter has already been fully disclosed as Step 1 for the calculation of the 
Fourier tr£insform, and the wavelet inverse filter or reconstruction algorithm has likewise 
already been described as Step 2 for the calculation of the inverse Fourier transform. 
Application-specific quantization is done according to psychovisual or psychoacoustic 
thresholds. Quantization and storage are furthermore merged with wavelet filtering 
into a single binary cascade as before, and retrieval or transmission are likewise merged 
with reconstruction into another single binary cascade. The reconstruction algorithm 
is exact. Furthermore, source code for it may be transmitted along with compressed 
data since the coding is sufficiently abbreviated and low-level. 

As is well-known, useful compression by transform coding requires an efficient specifi- 
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cation of the basis elen .ts. This is especially advantageous for new wavelets: a 
top arrow of generation g labeled by the two-by-two integral matrix A is specified by 
g bits, namely, by writing A uniquely as the matrix product of factors U or T. For 

example, the matrix = ^ 7^ factors as 

(lo V) = {1 ?)(o 1)' 

so the coefficient for the wavelet is coded by the binary sequence TUUUTTU, 

The method is also especially well-suited to progressive picture build-up or other it- 
erative refinement as follows. The Farey quadrature determines a linear ordering 
on the set Q of quadrature points in the circle C as before. The ordered subset 
2,±1} Q Q Q C is put into bijective correspondence with the collection of all 
arithmetic arrows (where \ denotes the set- theoretic relative complement); as illustrated 

in Figures la-lc, the arithmetic arrow labeled by ^ = corresponds to the com- 

plex number 

[ttZZZl ^ for a top arrow; 
gEgrllEfj e C, for a bottom arrow; 
i G Cy for the doe. 

Thus, the ordering from the Farey quadrature on Q\{~i, ±1} determines a correspond- 
ing ordering on the set of all arithmetic wavelets. Storage, retrieval, transmission, and 
reconstruction of wavelet coefficients is accomplished at a specified input or output data 
resolution in this canonical ordering since energy compacts to the lesser wavelet coef- 
ficients. In the case of real-time compression and transmission or other manipulation 
not requiring retrieval, the method can be further improved by the manipulation of 
previously computed wavelet coefficients in parallel with the calculation of subsequent 
ones. 



Complex and Multi-Dimensional Data 

There is nothing remarkable about the extension of the invention to complex-valued 
input data f(0) = u{d) -h i v{9), where u and v are real-valued functions on the circle 
C Concentrating for definiteness on arithmetic wavelets, one simple approach is to 
separately and independently transform u and 

-Yi^^ ^^^^^ ^(^) ^ X^/>i ^aW, 

using the techniques already described in order to calculate 

Another direct approach is to allow complex values for the input function / and solve 
for complex wavelet coefficients using the same algebraic formulas as before. 
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rurning now to multi-. lensional input/output data, suppose th there are R > 1 
input functions each of which has M > 1 independent variables, and define the torus 
C^^ to be the M-fold Cartesian product 

= C X . . . X C . 

M times 

Coordinates on C^'^ are given by M-tuples ^= (^i, • • . , 9m) of angles, and the input 
data is specified by an i?-tuple of mult i variable real- or complex- valued functions 

Using the bijection described in the previous section on data compression, each quadra- 
ture point Z e 2, ±1} uniquely determines a corresponding arithmetic wavelet 
'^zi^)] by convention, one also defines i9-t(0) = 'd±i{0) — 1. Concentrating on arith- 
metic wavelets for definiteness, define the corresponding multiwavelet 

As before, the input function F must be normalized. To this end, define a trigonometric 
monomial to be a product of the form X^{ei)X^(62) - • • ^^(^m), where each X\9) is 
given by one of X^{e) = cos 6, X\e) = sin 0, or = 1. The values of F at 

the points in {— z, ±1}^ uniquely determine coefficients in a hneax combination u{§) 
of trigonometric monomials so that u{e) = F{9) for each point of Just 
as in the one-dimensional^ wavelet filtei^ already discussed, now one must consider the 
normalization F(9) = F(9) - i/(0), so F is zero at each point of {-i, dbl}^. 

Adopting the Farey quadrature in each factor circle C of the torus C^, there is an 
induced Jexicographic ordering on the set Q^^\{~i,±l}^ C of potential sample 
points Z e C^^ , One separately approximates each coordinate function 

f{9) « i9^((9), for 2 = 1, 2, . . . , iJ, 

of F{9) as a finite linear combination of multiwavelets with application-specific zonal 
sampling for optimum energy compaction; this calculation may employ least-squares or 
other regularization in each factor as well as renormalization as for the one-dimensional 
wavelet filter. The methods are again elegantly implemented in practice as M nested 
binary cascades. The final processing for the cascades of arrow-structures in the cal- 
culation of multi-dimensional Fourier transforms may involve adding suitable linear 
combinations of trigonometric monomials based on further sampling, or one may sim- 
ply truncate with no final processing at all, as in the one-dimensional case. 

These same remarks apply verbatim to fan and modular wavelets. 



Mathematical Basis 

"The decorated Teichmliller space of punctured surfaces" , Communications in Mathe- 
matical Physics 13, pp. 299-339 (1987), written by the author of this patent applica- 
tion, began the study of certain abstract geometric coordinates called lambda lengths 
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in the context of hype Aic geometry on Riemann surfaces. Th^ leas developed in 
this and other related papers have found applications in high energy physics. A more 
recent publication in this series of about 20 papers is entitled "Universal constructions 
in Teichmiiller Theory", Advances in Mathematics 98, pp. 143-215 (1993). Here it 
was shown how to generalize lambda lengths to provide coordinates on a certain space 
of homeomorphisms (i.e. continuous bijections whose inverses are also continuous) of 
the circle, where there is one generalized lambda length for each arithmetic arrow. 
In further recent work described in "The Lie algebra of homeomorphisms of the cir- 
cle", Advances in Mathematics 140, pp. 282-322 (1998), written jointly with Feodor 
Malikov, there is described a representation of the coordinate deformations of these 
generalized lambda lengths as explicit vector fields on the circle. These vector fields 
are called "normalized elementary vector fields", and their study led to other vector 
fields on the circle called "fans" and "hyperfans". These three families of vector fields 
correspond, respectively, to the arithmetic, fan, and modular wavelets described here 
via the identification of the vector field f{0)§^ with the function f{9)\ there was, how- 
ever, no discussion of wavelets, sampling of data, approximation of functions, or any 
algorithms in the published literature. 

Each arithmetic wavelet is once-continuously difFerentiable on the circle, compactly 
supported, and localized in space. Arithmetic wavelets are not compactly supported 
in frequency, but instead, the frequency profile of an arithmetic wavelet is given alge- 
braically by the formula for used in Step 2 in the calculation of Fourier transforms; 
a non-compactly supported localization in frequency follows directly from this. The 
asserted formula for can be derived without much diflficulty but in several cases 
directly from Figures 3a-3e integrating twice by parts the usual expression for Fourier 
coefficients using the fact 'Oa is once-continuously differentiable. More difficult to derive 
is the formula given for the exponential functions e*"^ in terms of arithmetic wavelets 
which was used in Step 1 of the calculation of inverse Fourier transforms; this subtle 
computation with lambda lengths is given in "The Lie algebra of homeomorphisms of 
the circle" . 

Similar remarks apply to fan and modular wavelets, but fan wavelets are only contin- 
uous (they are not differentiable), and modular wavelets are not even continuous. 

In order to explain the failure of orthonormality of arithmetic wavelets and illuminate 
the regularization used in the arithmetic wavelet filter, fix some matrix A labeling 
an arithmetic arrow whose underlying chord has initial point q in the circle C. The 
sequence V^A of matrices label the arithmetic arrows whose underlying chords also 
have q as initial point, where n denotes an arbitrary integer. There is one infinite 
linear dependence 



for each q G Q, and this explains the additive degree of freedom in the calculation of 
euA and erA in the wavelet filter which is handled by regularization as was discussed 
before. This additive degree of freedom is also manifest in the initial specification 
6/ = 0 of wavelet coefficient on the doe in the wavelet filter. 

In order to better understand the Farey quadrature, it is convenient to transform the 



0 = 5^7?c/n^. 
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disk D to the upper h plane U = {u + iv : v > 0} hy means of '' transform 



This function if is a homeomorphism on the interior of D which maps the unit circle 
C bijectively to the real-axis {u iv : v = Q) together with an additional point oo 
at infinity, where K{-1) = 0, K{-i) = 1, and K{-hl) = oo. Furthermore, as ^ € Q 
varies over all quadrature points, the real number K{q) varies over all rational numbers, 
where the rational number ^ corresponds to the quadrature point e Q, In order 
to illustrate the Farey tesselation itself in U, whenever Qi,q2 G Q are the endpoints 
of the chord underlying an arithmetic arrow in Z>, draw the semi-circle in U with real 
endpoints ^1,^2; in the degenerate case that one of the quadrature points is infinite, 
say the point q2 ~ 00, then draw the vertical ray in U with endpoint qi . Applying this 
procedure to each arithmetic arrow produces the classical model of the Farey tesselation 
in U that is indicated in Figure 11. This figure illustrates the relationship between 
the Farey tesselation and the indicated circle packing, where the circle in the figure 
that is tangent to the real-axis at the rational point p/q has diameter q~^. For further 
information on number- theoretic aspects of this classical figure, see, for instance, Ford's 
book "Automorphic Functions", Chelsea (1972). 



In order to adequately disclose the methods, it is the purpose of this paragraph to 
include source codes for their implementations in the computer language C. Source 
code is presented for each of the following two implementations. 

o awft is the implementation of a preferred embodiment of the method 
for computing the Fourier transform of real-valued input data defined on 
the circle. The algorithm depends upon a bandwidth N, tolerance EPS, 
generation cut-off NVAN, and minimum generation MING as described 
before. A further flag SMODE=l determines that the output data will 
include the 0,+!,-! Fourier coefficients, and SMODE=0 otherwise. (The 
only configuration in either code which requires any library routines is 
SMODE=l, which requires math.h.) 

o awift is the implementation of a preferred embodiment of the method 
for computing the inverse Fourier transform of specified complex- valued 
Fourier coefficients defined in a specified bandwidth N. The output data 
is controlled by a single parameter SCAT, which is related to the param- 
eter SCALE discussed before by SCAT=[exp sinh " (-1)(1/SCALE)] ' 2. 
(SCAT varies inversely with SCALE and may be compared with conve- 
niently derived approximate quantities.) 

The comparison of the subroutines tgener in awft and in awift illustrates the easy 
transition between source codes for real and for complex data. The corresponding 
complex-awft and real-awift are therefore easily derived from the included source codes. 



K :D-^U 



z ^ % 



.2 + 1 



Source Codes 
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Furthermore, the met. J for compression of one-dimensional in] data amounts to 
the first step of awft, then quantization/de-quantization, then the second step of awift. 
The included source codes together thus also implement the method for compression 
of one-dimensional data since source code for it may be extracted from the included 
source codes. Finally, it is straight-forward to write the further driving routine required 
for multi-dimensional data compression, so the source codes included with this patent 
apphcation are sufficient to readily implement the method for multi-dimensional data 
compression as well. 

The included implementations have not been optimized. (For instance, straight-forward 
optimization has improved run-time by as much as 70% over the included awft code; 
on the other hand, optimization destroys intelligibility of the coding in relation to 
the underlying algorithm, and this explains the inclusion of unoptimized source code.) 
Another substantial improvement in run-time over the included source codes can be 
achieved by factoring so as to replace various manipulations of complex numbers with 
manipulations of integers as discussed before. 

Moreover, for clarity, the included source codes do not take advantage of renormal- 
ization, but for completeness, to each of awft and awift is appended source code to 
replace the subroutines tgener and bgener by corresponding subroutines which do em- 
ploy renormalization. 
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//This is an implementaiion in C of a preferred embodiment 
//of the invention for the calculation of the Fourier transform 



#include <stdio.h> 
#include <math.h> 
#define GEN 1000 
#derme N 50 
#define COR 1 
#define EPS l.e-12 



//math.h is required only for the calculation of 0,+l,-l Fourier coeffs 
//the maximum allowed generation 
//output bandwidth 
//irrelevant parameter 

//tolerance, where contributions to Fourier coeffs of norm less than 
//EPS are regarded as negligible 



^define PI 3.141592653589793 

#define NVAN 1 //terminate the recursion when NVAN consecutive 

//generations have contributed negligibly 

#define MING 15 //calculate at least MING generations 



#defme SMODE 1 



//SM0DE=1 includes the calculation of 0.+l,-l Fourier coeffs 
//SMODE=0 is faster and does not require math.h 



struct complex { //a complex number 

double x; 
double y; 

); 



struct trip{ 

double alp; 
double bet; 
double gam; 

); 



//convenient to have a real triple 



struct edge{ 

int a,b,c,d; 

double e,aip.bet,gam; 



//an arrow-structure 



int count=0; 



//a running count of the total number of samples of input data 



struct complex fourier[2*N+l]; //the complex Fourier coeffs stored in the order 

//n=-N,...,.l Al N. so fourier(n) is (N+n)th coeff 
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double abar,bbar,cbar; //coeffs in trig function abar*cos theta + bbar*sin theia + cbar 

//which agrees with input function for theta = -PI, -PI/2, 0 

void njain(void){ 

void genertopCstruct edge *.int.int,double); //recursive routine for top of circle 
void generbot(struct edge *,int,int4oubie); //recursive routine for bottom of circle 

void normalout(void); //output normalization routine 

void nonnalin(void); //input normalization routine 

double fbar(int.int); //fbar(p.q) returns f(theta) - (abar*cos theta + bbar*sin theta + cbar) 
//where theta is the argument of the complex number (p-iq)/(p+iq) 

struct edge doe[l]; //initial edge for recursive calls 

inti; 

normalinQ; //this calculates the normalization parameters abar, bbar, cbar 

doe->a=:doe->d=l ; //initialize the matrix of the doe to the identity 

doe->b=doe->c=0; 

doe->alp=doe->bet=doe->gam=doe->e=0; //initialize wavelet and lagged trig coeffs 

genertopCdoe. 1 .0.fbar(- 1. 1)); //begin recursion on top of the circle 

generbot(doe, 1 .0,0.); //begin recursion on bottom of the circle 

normaloutO; //this divides each output Fourier coeff by an overall factor PI 

//which was suppressed during the recursion and determines 
//negative from positive Fourier coeffs using reality of input 

printf("with EPS=%e, NVAN=%d, COR=%d, and MING=%d. the count is: %d\n", 
EPS.NVAN.COR,MING,count); 

for(i=0;i<=:N;i-M-) 

printf("c%3d = %e-f(i)%e\n",i,fourier[N+i].x.fourier[N+i].y); 



) 



void genertop(struct edge *p,int g jni envy,double fc){ //recursion for top of circle 

//input pointer to a current edge p of generation g, where there 
//have so far been envy consecutive generations of negligible 
//contributions, and fc is the current normalized input sample 
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struct edge twofer[2]; 



// the descendant edges of p 



struct complex temp; 



void genertop(struct edge *,int,int,double); //recursive function 



struct complex tgener(struct edge *,struct edge *,doub!e); 



//tgener(p.twofer,0 stuffs iwofer with the descendants of p 
//using the updated normahzed input sample f returning the 
//next two normalized samples as its real and imaginary pans 



ini charles(struct edge *); //charles(p) updates the Fourier coeffs with the 

//contribunon from edge p returning 1 if they 
//are negligible, and returning 0 otherwise 

void tprune(struct edge *); //this implements the processing of a terminal edge 

if(envy>=NVAN && g>= MING)( //if appropriate, then terminate cascade 
tprune(p); 
return; 
) 

temp=tgener(p,twofer,fc); //otherwise, generate descendants of p 

if(charles(twofer)*charles(twofer+l)) //if both descendants contribute negligibly, 

envy++; //then increment envy 



else 



envy=0; 



//otherwise reset envy to zero 



if(g>=GEN){ 



//problem: non-convergence after GEN generations 



printfCfall through with gen=%d for a.b,c,d=%d,%d,%d,%d and e=%e\n". 
g,p->a,p->b,p->c,p->d,p->e); 

return; 



else( 



//otherwise, continue recursion with descendants 



genertop(twofer,g+ 1 ,envy,temp.x); 
genertop(twofer-f 1 ,g+ 1 ,envy,temp.y); 



void generbot(struct edge *p,int g,int envy,double fc)( 



//recursion for bottom of circle 
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struct edge twofer[2]; 



//this is entirely analogous to genertop 



struct complex temp; 

void generbot(struct edge *.int,int.doubIe); 

struct complex bgener(struci edge *,struct edge *,double); 

int charles(struct edge *); 

void bprune(struct edge *); 

if(envy>=NVAN && g >= MING)( 
bprune(p); 



temp=bgener(p,twofer,fc) ; 

if(charles(twofer)*charles(twofer+i)) 
envy++; 

else 

envy=0; 
if(g>=GEN){ 

printfC'fall through with gen=%d for a,b.c.d=%d,%d,%d.%d and e=%e\n*', 
g»p->a,p->b,p->c,p->d.p->e); 

return; 
) 

elsef 

generbot(twofer,g+ 1 ,en vy,temp.x); 
generbot(twofer-t-l ,g+ 1 .envy,temp.y); 
) 



struct complex makecpx(intjnt); //input p,q returns the complex number (p-iq)/(p+iq) 
struct complex mult(struct complex^struct complex); //complex multiplication 



return; 



} 



int charles(siruci edge ♦ p){ 



//charles(p) updates the Fourier coeffs using p 
//and returns a control parameter flag=l if 
//contribution is negligible and flag=0 otherwise 



int n,flag; 



int a,b,c,d; 
double e; 



//local variables for matrix values 
//local variable for wavelet coeff 



struct complex pzl,p2r,pzf,p2b; 



//l,r»f,b labels !eft,right,froni,back vertices 
//of Farey tessetation near the edge p. and 
//in a loop below on index n, p2x=(pxm)^n 



struct complex plm,prm,pfm,pbm; 
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//for x=l,r,f,b 



struct complex tmp,mul,cn; 



int sixupc(ini,int,int»int,doubIe); 



//this updates the Fourier coeffs 0,+l,-l 
//provided SMODE=l 



a=p->a; 
b=p->b; 
c=p->c; 
d=p->d; 
e=p->e; 



//stuff local values 



p2r.x=b-d; 



//set-up for loop over Fourier coeffs 



p2r,y=c-a; 

prm=makecpx(b-d,a-c); 
pzr=mult(p2r,pzr); 

p2l.x=b+d; 
p2l.y=-c-a; 

plm=makecpx(b+d,a+c); 
pzi=mult(p2l,p2l); 



pzf.x=d; 
pzf.y=-c; 

pfm=makecpx(d,c); 
pzf=inuU(p2f,p2f); 



p2b.x=b; 
P2b,y=-a; 

pbm=makecpx(b,a); 
p2b=mult(p2b,p2b); 

flag=l; //initialize flag 

tmp,x=0; 

if(SMODE) 

flag=sixupc(a,b,c.d.e); //update 0.+ 1 1 modes and re-set flag to 0 if any 
//of these contributions is non-negligible 

for(n=2;n<=N;n-M-) { //loop over positive Fourier coeffs n> 1 

tmp.y=e/(n-n*n*n); 

p2r=mult(pzr,piTn); //recursively calculate powers 

P2l=mult(p2l.pim); 
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pzf=:mult(i.. ..pfm); 
pzb=muit(pzb,pbm); 



cn.x=-pzr.x+2*pzf.x+2*p2b.x-p2l.x; 
cn.y=-pzr.y+2*pzf.y+2*p2b.y-pzi.y; 

mul=mult(cn.tmp); 

fourier(n+N].x+=mul.x; //update nth Fourier coeff 
fourier[n+N].y+=muI.y; 

if(nag==:l && mul.x*mul.x+niuLy*muLy>=EPS) //if contribution to nth Fourier 
^•^5=0; //coeff is non-negligible. 



//then set flag to zero 



retum(flag); 



int sixupc(int a,int b.int c.int d,double e){ 



//this updates the 0,+ !,- 1 Fourier coeffs 
//using the input matrix entries a,b,c,d 
//and wavelet coeff e 



double tp.tm,tr,tl; //p,ni,r,l labels the tenminaUinitiaI,rightJeft vertices of Farey 

//tesselation near edge, and tx denotes corresponding angle 

double dtemp; 



double tb; 



//tb=l on the top and tb=-l on the bottom of the circle 



int fl; 



//this is the local version of the flag in Charles and is die 
//value returned to Charles by sixupc 



struct complex makecpx(int,int), temp; 



fl=l; 



//set the flag to 1 



tb=:l,; 
if(b+c<0) 



//default to the top of the circle 
//if on the bottom of the circle 
//then re-set tb=-l. 



tb=-l.; 



//calculate the various angles 



temp=makecpx(d.-c); 
tp=acos(tb*temp.x); 



temp=makecpx(b,-a); 
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tm=acos(tb* temp, x^, 

temp=makecpx(b+d,-(a+c)); 
tl=acos(tb*temp.x); 

temp=makecpx(b-d,c-a); 
tr=acos(tb*temp,x); 

dtemp=e*( //contribution to 0th Fourier coeff 

2.*tp*(c*c+d*d)+2.*tm*(a*a+b*b) 
-tl*((a+c)*(a+c)+(b+d)*(b+d)) 
-tr*((a-c)*(a-c)+(b-d)*(b-d)) 

); 

if(dtemp*dtemp>=:EPS) //if this contribution is non-neghgible, 
fl=0; //then re-set the flag 

fourier[N].x+=dtemp; //up-date the 0th Fourier coeff 

temp.x=e*( //contribution to real part of 1st 

tp*(c*c-d*d)+tm*(a*a-b*b) //Fourier coeff 

+tr*((b-d)*(b.d).(a-c)*(a-c))/2, 
+tl*((b+d)*(b+d)-(a+c)*(a-K:))/2, 

); 

fourier[N+l].x+=temp.x; //update real part of 1st Fourier coeff 



iemp.y=e*( //contribution to imag part of 1st 

2 *(tp*c*d+tm*a*b) //Fourier coeff 

-tr*(b-d)*(a.c) 
-tl*(b-Ki)*(a+c) 

); 

//update imag part of 1st Fourier coeff 

//if this contribution is non-negligible, 
//then re-set the flag to zero 



fouriertN+l).y+=temp.y; 

if(temp.x*temp.x+temp.y*temp.y>=EPS) 
n=0; 

retum(fl); 



struct complex muU(struct complex u,struct complex v){ //complex multiplication 

struct complex temp; 
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temp.x=u.x*v.x-u.> .y; 

temp.y=u.x*v.y+v.x*u.y; 



retum(temp); 



struct complex makecpx(int p,int q){ 



//input p.q returns the complex number (p+iq)/(p+iq) 



struct complex temp; 
double den; 

den=p*p+q*q; 

temp.x=(p*p-q*q)/den; 
temp.y=:-2*p*q/den; 

retum(temp); 



double fz=0..fi=0.,fo=0.; 
double alp,bet,gam; 
int n; 

for(n=0;n<=N;n-M-) { 

fourier[N+n],x=fourier[N+n].x/PI; //include factor PI suppressed before 
fourier[N+n]-y=fourierrN+n].y/PI; 

fourier[N-nl,x=fourier[N+n].x; //calculate negative from positve 

fourierfN-n].y=-fourier[N+n].y; //Fourier coeffs 

) 



void normalout(void){ 



//output nomialization 



void normaiin(void){ 



//input normalization 



double fz=0..ri=0.,fo=0.; 



//values of input function for angles -PI, -PI/2, 0 



int i; 



double f(int,int); 



//f(p»q) returns the un-normalized input sample 
//at the angle corresponding to the point on the 
//circle given by the complex number (p-iq)/(p+iq) 



for(i=0;i<=2*N;i-H-) 

fourier[i].y=fourier[i].x=0; 
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fz=f(0,l); 
fi=f(l,0); 
fo=f(l,l); 

abar=(ri-fz)/2; //abar*cos theta + bbar*sin theta -h cbar takes the same values 

cbar=Cfi+fz)/2; //as f for the angles theta=-PI, -PI/2, 0 

bbar=cbar-fo; 



fourier[N].x=:cbar*PI; 
fourier[N].y=0.; 

fourier[N+ 1 ] .x=abar*PI/2; 
fourier[N+ 1 ].y=-bbar*PI/2; 



//stuff the 0th Fourier coeff scaled by PI to 
//recover it after re-scaling by 1/PI in normalout 

//stuff the Isi Fourier coeff scaled by PI to 
//recover it after re-scaling by 1/PI in normalout 



double fbar(int p, int q){ //input p.q to fbar produces the nomalized value 

struct complex temp; //of the input data at the point (p-iq)/(p+iq) 

struct complex makecpx(int»int); 
double f (intent); 
temp=makecpx(p.q); 

retumCf(p,q) - (abar*temp.x + bbar*temp.y-i-cbar)); 

) 

double f{int p.int q) { //f is the input data, where f(p,q) is the 

double fsinm(int.int.int); //value taken at (p-iq)/(p+iq) 

retum(fsinm(p,q.6)); 

} 

double fsinm(int p,int q, int m){ //for example, the input function is 

ii^t /^e^e taken to be sin(6*theta) 

struct complex zeta,meta; 

struct complex makecpx(int,int),mult(struct complex .struct complex); 

zeta=makecpx(p,q); 

meta.x=l.; 

meta.y=sO; 

for(n=l ;n<=m;n++) 

meta=mult(2eta»meta); 

retum(meta.y); 

) 



struct complex tgener(struct edge *pc,struct edge *pn,doubIe fc){ 

//tgener(p. twofer,fc) stuffs twofer with the descendants 
//of p for the top of the circle, where fc is the current 
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//normalized input sample, and tgener returns a complex nuiuoer. 
//the real and imag pans being the normalized input values 
//required for the least-squares regularization 



double fbar(int,int); 

int a.b»c.d,bpd,apc; 

double e,alp,bei,gam; 

double axJ,ayl.ax2,ay2; 
double bx 1 ,by 1 ,bx2,by2 ; 
double cxl,cyl,cx2,cy2; 

double Sigl.taul;sig2,tau2,chi; 

double eu.et; 

double fn.fnp; 



//fbar(p,q) returns the normalized input values 

//local variables for matrix entries, bpd=b+d. apc=:a+c 

//local variables for wavelet and lagged trig coeffs 

//index 1,2 refers to first.second (counter-clockwise) 
//descendants* prefix a.b.c refers to alp.bet,gam 
//update procedures below define indices x,y 

//quantities used in the least-squares regularization 

//wavelet coeffs of first and second descendant edges 

//normalized input samples 



double vlc,vlu,vlt.v2c,v2u.v2i; 

//coeffs in the inhomogeneous linear expression vic+viu*eu+vit*et, i=l,2. 
//of ongoing approximation to normalized input at the respective sample points 
//((2b+d)+i(2a+c))/((2b+d)-i(2a+c)), ((b+2d)+i(a+2c))/((b+2d)-i(a-i-2c)) 



struct complex temp; 
struct edge *pnp; 



pnp=pn; 
pnp-Hf-; 

a=pc->a; 
b=pc->b; 
c=pc->c; 
d=pc->d; 
apc=a+c; 
bpd=b^-d; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 



//initialize pointers 



//initialize matrix entries 



//stuff matrix of first descendant 



//stuff matrix of second descendant 
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alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



//initialize lagged trig coeffs 



e=pc->e; 



//initialize wavelet coeff 



temp.x=fn=:fbar(-(b+bpd)»a+apc); //store new normalized input samples as the 
temp.y=fnp=fbar(-(d+bpd),c+apc); //real and imag parts of temp to return 



axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 
ayl=2*(b*b-a*a); 

bxl=bet+(c*d-a*d-b*c-a*b)*4*e; 
byl=4*a*b; 

cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; 
cyl=-2*(a*a+b*b); 

ax2=alp+4*(d*d-c*c)*e; 
ay2=2*(c*c-d*d); 

bx2=bet+8*e*c*d; 
by2=-4*c*d; 

cx2=gam-4*e*(c*c-Ki*d); 
cy2=2*(c*c-Kl*d); 

chi=2*e-K //chi=eu-et 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

) 

//prepare for stuffing alp,bet,gam in case c>=a 



axl=alp+4*e*(a*a-b*b); 
ayl=2*(b*b.a*a); 

bxl=bet-8*e*a*b; 
byl=4*a*b: 

cxl=gam+4*e*(a*a+b*b); 



count+=2; 



//update count of total sample points 



if(a>c){ 



//prepare for stuffing alp.bet»gam in case a>c 
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cyl=-2*(a-a+b*b); 

ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c+c*d)*4*e; 
by2=-4*c*d; 

cx2=gam-Ka*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; 
cy2=2*(c*c+d*d); 



chi=-2*e+( //chi=eu-et 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gani) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

} 

//in any case, calculate the coeffs vic,viu,vit, for i=l,2 
vlc=:-fn+cxl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl 
+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

vlt=cyl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl 
+2*byl*(b+bpd)*(a+apc))/((b+bpd)*(b4.bpd)+(a+apc)*(a+apc)); 

V 1 u=8,/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

v2c=-fnp+cx2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2 
+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2 
+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

sigl=vlc+chi*vlu; 
taul=vlt+vlu; 

Sig2=v2c+chi*v2u; 
tau2=v2i+v2u; 

et=-(sig I ♦tau 1 +sig2*tau2)/(tau 1 *tau 1 +tau2*tau2); //calculate et and eu 
eu=chi+et; 
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pn->e=eu; 
pnp->e=el; 



//stuff new wavelet coeffs 



pn->alp=ax 1 +ay 1 *et; 
pn->bet=bxl+byl*et; 
pn->gam=cx I +cy 1 *et; 



//stuff first lagged trig coeffs 



pnp->alp=ax2+ay2*eu; 



//stuff second lagged trig coeffs 



pnp->bet=bx2+by2*eu; 
pnp->gam=cx2+cy2*eu; 

return (temp); 

) 



struct complex bgener(struct edge *pc,struct edge *pn.double fc){ 



/^gener(p,twofer,fc) stuffs twofer with the descendants 
//of p for the bottom of the circle, where fc is the current 
//normalized input sample and returns a complex number, 
//the real and imag parts being the normalized input values 
//required for the least-squares regularization 

//this is entirely analogous to tgener. and only the differences 
//will be noted here 



double fbar(int,int); 

int a,b,c,d,bpd,apc; 
double e,alp,bet,gam; 

double axi,ayl,ax2,ay2; 
double bxl.byl,bx2»by2; 
double cxl,cyl,cx2.cy2; 

double sigl,taul,sig2,tau2.chi; 
double vlc.vlu.vlt.v2c,v2u»v2t; 

double eu.et,fn,fnp; 

struct complex temp; 

struct edge ♦pnp; 

pnp=pn; 
pnp++; 
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a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 

apc=a+c; 
bpd=b+d; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 

pnp->a=d; 
pnp->b--c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc*>alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 

remp.x=fn=fbar(a+apc,b+bpd); //differs from tgener in that 
temp.y=fnp=fbar(c+apc»d+bpd); //(-q,p) replaces (p.q) 

count+=2; 

if(a>c)( 

axl=alp-Kd*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 
ayl=:2*(b*b-a*a); 

bxl=bet+(c*d-a*d-b*c-a*b)M*e; 
byl=4*a*b; 

cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; 
cyl=-2*(a*a+b*b); 

ax2=alp+4*(d*d-c*c)*e; 
ay2=2*(c*c-d*d); 

bx2=bei+8*e*c*d; 
by2=-4*c*d; 

cx2=:gam-4*e*(c*c+d*d); 
cy2=:2*(c*c+d*d); 



//differs from tgener in that d,-c,-b,a replaces a,b,c,u 



//differs from tgener in that d,-c,-b,a replaces a,b,c.d 



//differs from tgener in that d»-c.-b,a replaces a.b,c,d 
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chi=2*e4-( 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 



ax 1 =alp+4*e*(a*a-b*b); 
ayl=2*(b*b-a*a); 

bxl=bei-8*e*a*b; 
by!=4*a*b; 

cxl=gam+4*e*(a*a-»-b*b); 
cyl=-2*(a*a+b*b); 

ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c+c*d)M*e; 
by2=-4*c*d; 

cx2=gam+(a*a+b*b-2*(a*c-fb*d)-c*c-d*d)*2*e; 
cy2=2*(c*c+d*d); 



chi=-2*e+( 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

} 

vlc=-fn+cxl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl 
+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

vlt=cyl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ay 1 

+2*by I *(b+bpd)*(a+apc))/((b+bpd)*(b-fbpd)+(a+apc)*(a+apc)); 

vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 



else{ 
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v2c=-fnp+cx2+ 

(((d+bpd)*(d-*-bpd)-(c+apc)*(c+apc))*ax2 
+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2 
+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2t--8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

sigl=vlc+chi*vlu; 
taul=vlt+vlu; 

sig2=v2c-K:hi*v2u; 
tau2=v2t+v2u; 

et=-(sig I *tau I +sig2*tau2)/(tau 1 ♦tau 1 +Iau2*tau2) ; 
eu=chi+et; 

pn->e=eu; 
pnp->e=et; 

pn->alp=axl+ayl*et;; 
pn->bet=bx 1 +by 1 *et; 
pn->gam=cx 1 +cy 1 ♦et; 

pnp->alp=ax2+ay2*eu; 

pnp->bet=bx2+by2*eu; 
pnp->gani=cx2+cy2*eu; 

return (temp); 



tprune(stnict edge *p) { //this routine implements the processing of terminal 

//edge p in the cascade for the top of the circle 

double alp, bet» gam.e; //local copies of lagged trig coeffs and wavelet coeff 

double alpl.alp2,betl.bet2,gaml,gam2; //trig coeffs of the trigonometric function sigma at 

//the first and second (counter-clockwise) samples 

struct trip mob(siruci complex.struct complex, struct complex, 

doubIe,double,doubIe); 
//mob(zl,z2.z3,fl,f2,f3) returns the triple with entries alp,bet,gam 
//where aIp*cos iheta + bet*sin theta + gam takes the respective values 
// fl,f2,f3 at the points 2l,z2,z3 on the circle given as complex numbers 
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double fbarCint,int); 

void fixupCstruct trip,int»int,int,int,int); 

struct complex makecpx(int,int); 

void sixupp(struct tripjnt,int jntant^int); 



//returns normalized input values 

//adjusts the output Fourier coeffs using sigma-tau 

//input p,q returns (p-iq)/(p+iq) 

//adj usts the 0,+ 1 1 Fourier coeffs provided 
//SMODE=l 



int a.b,c,d; //local variables for the matrix entries 

struct trip tripl,trip2; //first and second trig triples 

a=p->a; //stuff local variables 

b=p->b; 

c=p->c; 

d=p->d; 

e=p->e; 

alp=p->alp; 

bet=p->bet; 

gam=:p->gam; 

if {a>c){ //final update of trig fields in case a>c 

alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); 

betl=bet+e*4*(c*(d-b)-a*(d+b)); 

gaml=gam+e*2*(a*(a+cHc*(a.c)+b*(b+d)-d*(d-b)); 

alp2=alp+e*4*(d-c)*(d+c); 
bei2=bet+e*8*c*d; 
gam2=gam-e*4*(c*c+d*d); 
) 

else{ //final update of trig fields in case c>=a 

alpl=alp+e*4*(a-b)*(a+b); 
betl=bet-e*8*a*b; 
gam 1 =gam+e*4*(a*a+b*b); 



alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); 
bet2=bet+e*4*(a*(d-b)+c*(b+d)); 
gam2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d)); 
) 

tripl=mob( //calculate sigma for first descendant 

makecpx(-(3*b+d)3*a+c). 
makecpx(-(2*b+d),2*a+c). 
makecpx(-(3*b+2*d).3*a+2*c). 
fbar(-(3*b+d),3*a+c), 
fi5ar(.(2*b+d),2*a+c), 
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fbarC-(3*b-h2*d).3*a+2*c) 

); 

trip2=mob( //calculate sigma for second descendant 

makecpx(-(2*b+3*d).2*a+3*c), 
makecpx(-(b+2*dXa+2*c), 
makecpx(-(b+3*d),a+3*c), 
fbar(-(2*b+3*dX2*a+3*c), 
fbar(-(b+2*d).a+2*c). 
fbar(-(b+3*d).a+3*c) 
); 

coun(+=6; //update count of total sample points 

tripl.alp=(tripl .alp-aipl)/COR; //difference of first trig functions 
trip 1 .bet=(trip 1 .bet-bet I )/COR; 
trip 1 .gam=:(trip 1 .gam-gam 1 )/COR; 

trip2.alp=(trip2.alp-a!p2)/COR; //difference of second trig funcdons 

trip2.bet=(trip2.bet-bet2)/COR; 

trip2.gam=(trip2.gam-gam2)/COR; 

rjxup(tripl,b+d,a+c.b.a,l); //adjust Fourier coeffs with first trig function difference 

fixup(trip2,d,c,b+d.a+c.l); //adjust Fourier coeffs with second trig function difference 

if(SMODE> { //if SMODE= 1 , then adjust 0,+l.-l Fourier coeffs as well 

sixupp(trip l,b-Ki.a+c,b.a. 1 ); 
sixupp(irip2,d,c,b+d,a+c. 1 ); 
) 



void bprune(struct edge *p){ //this routine implements the processing of terminal 

//edge p in the cascade for the bottom of the circle 

//this is entirely analogous to tgener, and only the 
//differences will be noted here 

double alp, bet, gam,e; 
double alp 1 ,alp2 ,bet 1 ,bet2,gam 1 ,gam2; 
struct trip mob(struct complex.struct complex, struct complex, 

double.double,double); 

double fbar(int jnt); 

void fixup(struct trip. intjnt4nt,ini.int); 

struct complex makecpx(int.int); 

void sixupp(struct trip,!nt.int.int,int.int); 

int a,b,c,d; 

struct trip tripl.trip2; 
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a=p->d; 



//differs from tprune in that d.-c,-b,a replaces a,b,c.d 



b=-(p->c); 

c=-(p->b); 

d=p->a; 

e=p->e; 

aip=p->alp; 

bet=p->bet; 

gam=p->gam; 

if (a>c){ 

alpl=aip+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); 
bet 1 =bet+e ♦A+Cc *(d-b)-a*(d+b)) ; 
gaml=gam+e*2*(a*(a+c)+c*(a-cHb*(b+d)-d*(d-b)); 

alp2=alp+e*4*(d-c)*(d+c); 
bet2=bet+e*8*c*d; 
gam2=rgam-e*4*(c*c+d*d); 
) 



alpl=alp+e*4*(a-b)*(a+b); 

betl=bet-e*8*a*b; 

gaml=gam+e*4*(a*a+b*b); 

alp2=aIp+e*2*(a*(a-c)+b*Cd-b)-c*(c+a)-Hl*(b+d)); 
bet2=bet+e*4*(a*(d-b)+c*(b+d)); 
gain2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d)); 
) 



tripi=mob( //arguments of mob differ from those in tprune 

makecpx(3*a+c,3*b+d). 
makecpx(2*a+c,2*b-Hl), 
makecpx(3*a+2*c,3*b+2*d), 
fbar(3*a+c,3*b+d). 
fbar(2*a+c,2*b+d), 
fbar(3*a+2*c,3*b+2*d) 

); 

trip2=mob( //arguments of mob differ from those in tprune 

makecpx(2*a+3*c,2*b+3*d), 
makecpx(a+2*c»b+2*d). 
makecpx(a+3*c.b+3*d). 
fbar(2*a+3*c.2*b+3*d), 
fbar(a+2*c,b+2*d), 
fbar(a+3*c,b+3*d) 

); 



else{ 
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tripl.alp=(mpl.alp+alpl)/COR; //differs from iprune in that alp,bet 

tripl.bet=(tripl.bet+betl)/COR; //is replaced by -alp.-bet 
trip 1 .gam=(trip 1 .gam-gam 1 )/COR; 

trip2.alp=(trip2.alp+aip2)/COR; //differs from tprune in that alp,bet 

trip2.bet=(trip2.bet+bet2)/COR; //is replaced by -a)p,-bet 
trip2.gam=(trip2.gam-gam2)/COR; 

fixup(tripl ,b+d,a+c.b.a,- 1); //last argument of fixup differs from tprune 

fixup(trip2,d.c»b+d.a+c,-l);Q //Jast argument of fixup differs from tprune 

if(SMODE)( //if SMODE=l, adjust the 0,+l,-l Fourier coeffs 

sixupp(tripl,b+d.a+c,b,a,-l); //last argument of sixupp differs from tprune 
sixuppCtrip2,d,c,b+d,a+c»-l); //last argument of sixupp differs from tprune 

) 



struct trip mob(struct complex z Instruct complex 22»struct complex 23. 

double fl,double f2.double f3){ 
//mob(zl,z2.z3.fl,f2,f3) returns the triple with entries alp.bet.gam 
//where aIp*cos iheia + bet*sin theta + gam takes the respective values 
// fl,f2,f3 at the points zl,z2,23 on the circle given as complex numbers 

double det.aa,bb,cc,dd; 
struct trip temp; 

aa=(zl.x-23.x); 
bb=(2l.y-z3.y); 
cc=(22.x-z3.x); 
dd=(z2.y-z3.y); 

det=aa*dd-bb*cc; //genera! principles guarantee that det is non-zero 

aa=aa/det; 
bb=bb/det; 
cc=cc/det; 
dd=dd/det; 

temp.aIp=dd*(fl-f3)-bb*(f2-0); 
temp.bet=aa*(f2-B)-cc*(fl-f3); 
temp.gam=f3-23.x*temp.alp-z3,y*temp.bet; 

retum(temp); 
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void fixupCstruct trip delt.int dm.int cm.int dp,int cp,int tb){ 

//fixup delt,dm.cm,dp,cp4b) adjusts the Fourier coeffs by adding the Fourier 
//coeffs of delt=T^C-T with support truncated to be the interval with endpoints 
//(dm+i *cm)/(dm-i *cm), (dp+i *dp)/(dp-i ♦cp) 

//let tm, tp be the corresponding angles 

// tb=+l on the top, and tb=-l on the bottom of the circle 

int n; 

double cpp,cmp,cpm»cmm; 
//in loop below on n. cpx=[cos(l+n)tx3/Cn+I), cmx=[cos(l-n)tx]/(l-n), for x=p,m 

double spm.smm. spp,smp; 
//in loop below on n, spx=[sin(l+n)tx]/(l+nX smx=rsin(l-n)tx]/(l-n). for x=p.m 

double crm.crp,smi»srp; 
//in loop below on n. cnn=[cos n*tx]/n. sm\= [sin n*tx]/n. for x=p,m 

double alp,bet.gam; //local variables for entries of delt 

double ctp, ctm, cnp,cnm; //in loop below on n. ctx=cos tx. cnx=cos n*tx, for x=p.m 

double stm,snin,stp,snp; //in loop below on n. stx=sin tx, snx=sin n*tx, for x=p,m 

double cbufp,sbufp,cbufm,sbufni; 
struct complex makecpx(ini,int),temp; 



alp=delt.alp; //stuff local variables 

bci=delt,bet; 

gam=dcU.gam; 



temp=makecpx(dm,-cm); 

ctm=tb*temp.x; 

stm=tb*temp.y; 

temp=makecpx(dp.-cp) ; 

ctp=tb*temp.x; 

stp=tb*temp.y; 

cnp=:ctp*ctp-stp*stp; 
snp=2*stp*ctp; 

cnm=ctm*ctm-stm*stm; 
snm=2*stm*ctm; 
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for(n=2;n<=N;n++) { //loop over Fourier coeffs 

cbufp=ctp*cnp-stp*snp; //update using formulas for cos 
sbufp=stp*cnp+ctp*snp; //or sin of a sum of angles 
cbufm=ctm*cnm-stm*snm; 
sbufm=stm*cnm+ctm*snm; 

cpp=:cbufp/( l+n); //update for p 

cmp=(cbufp+2*stp*snp)/( 1 -n); 
spp=sbufp/(l+n); 
smp=(sbufp-2*ctp*snp)/(l -n); 

cpm=cbufm/( l+n); //update for m 

cniin=(cbufm+2*stm*snm)/(l-n); 
spm=sbufn[i/( l+n); 
smm=(sbufm-2*ctm*snm)/( 1 -n) ; 

crp=2.*cnp/n; 
crm=:2.*cnm/n; 
srp=2.*snp/n; 
snn=2.*snm/n; 

cnp=cbufp; //update for next iteration 

snp=sbufp; 

cnm=cbufm; 

snm=sbufm; 

fourier[N+n].x+=( //adjust real part of nth Fourier coeff 

+alp*(smp+spp-smm-spm) 

+bet*(-cpp-cmp+cpm+cmm) 

+gani*(sip-srm) 

)/4.; 

fourier[N+n].y+=( //adjust imag part of nth Fourier coeff 

+alp*(cpp-cmp-cpm+cnim) 

+bet*(spp-smp-spm+smm) 

+gam*(crp-crm) 

)/4.; 



return; 

) 



void sixupp(struct trip triple,int dm^int cm,int dp,int cpjnt tb)( 

//sixup(delt,dm,cm.dp,cp,tb) adjusts the 0,+ l.-l Fourier coeffs by adding the 0,+ l,-l Fourier 
//coeffs of deltsT^C-T with support truncated to be the interval with endpoints 
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//(dm+i *cm)/(dm- i *cm). (dp+i -^dpVCdp-i *cp) 

//let im, tp be the corresponding angles 

// tb=+l on the top, and ib=-l on the bottom of the circle 

double alp,bet,gam; //local variables for entries of delt 

double com,sini»cop,sip,tm,tp; //cox=cos tx, six=sin tx, for x=p.m 

struct complex temp; 

struct complex makecpx(int,int); 

temp=makecpx(dm,-cm); //calculate com.sim.tm 

com=tb*temp-x; 

sim=tb*temp.y; 

tm=acos(com); 

temp=makecpx(dp,-cp); //calculate cop.sipap 

cop=tb*temp-x; 

sip=tb*temp.y; 

tp=acos(cop); 

alp=tb*triple.alp; //adjust alp.bet for the bottom of the circle 

bet=tb*triple.bet; 

gam=triple.gam; 

fourier[N].x+=( //update 0th Fourier coeff 

gam*(tp-tm)+a!p*(sip-sim)-bet*(cop-com) 
)/2.; 

fourier[N+l J.x+=( //update real pan of Isi Fourier coeff 

alp*(tp-tm)/2. 
+gam*(sip-sim) 

-bei*(cop*cop-sip*sip-com*com+sim*sim)/4. 

+a]p*(cop*sip-com*sim)/2. 

)/2.; 

fourier[N+l].y+-( //update imag part of 1st Fourier coeff 

-bet*(tp-tm)/2, 
+gam*(cop-com) 

+alp*(cop*cop-sip*sip-com*com+sim*sim)/4. 

+bet*(cop*sip-com*sim)/2. 

)/2.; 

) 
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//Here are the subroutines replacing igener and bgener above in an 
//implementation which takes advantage of re normalization 



struct complex igener(struct edge *pc,struct edge *pn4oubIe fc) 
{ 



double fbar(int.int); 
int a,b.c»d,bpd.apc; 
double e,alp,bet,gam; 
double axl.ayl,ax2,ay2; 
double bxl»byl.bx2,by2; 
double cxl,cyUcx2,cy2; 
double sigl,taul.sig2,tau2,chi; 
double vlc,vlu,vU,v2c.v2u»v2t; 
double deriv(int,int4nt,int,int,int); 
double dc,dn,dnp; 
double eu,et,fn,fnp; 
struct complex temp; 
struct edge *pnp; 
extern int count; 

pnp=pn; 
pnp++; 



a=pc->a; 
b=pc->b; 
c=pc->c; 
d:=pc->d; 

apc=a+c; 
bpd=b+d; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 



>> 



pnp->b=bpd; 

pnp->c=c; 

pnp->d=d; 
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a!p=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 



temp.x=:fn=fbar(-(b+bpd).a+apc); 
temp. y=fnp=fbar(-Cd+bpd).c+apc); 
temp,n-0; 

count+=2; 

dc=deri v(a, b,c,d,- 1,1 ); 
dn=deri v(a,b,c,d,- 1,2); 
dnp=deriv(a»b,c,d.-2,l); 

fc=fc/dc; 
fn-fn/dn; 
fnp=fnp/dnp; 



if(a>c){ 

axl=alp+4.*e; 

bxl=bet-4.*e; 

cxl=gam; 

ax2=alp+4.*e; 

bx2=bet; 

cx2=gam-4.*e; 

chi=2.*((fc-gam-bet)/4.+e); 

} 



else( 



axl=alp+4.*e; 

bxl=bet; 

cxl=gam-h4.*e; 

ax2=salp+4.*e; 

bx2=bet+4.*e; 

cx2=:gam; 

chi-2.*((fc-gam-bet)/4.-e); 
) 

V 1 c=-fn+cx 1 +(4, *bx 1 -3.*ax 1 )/5. ; 

vU=-4./5.; 

vlu=8./5.; 
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v2c=-fnp+cx 2+(4. * o.^ 2+3 *ax2)/5 . ; 

v2u=4./5.; 

v2t=-875.; 

sigl=vlc+chi*vlu; 
taul=vlt+viu; 

sig2=v2c+chi*v2u; 
tau2=:v2t+v2u; 

et=-(sig 1 *tau 1 *( 1 )+sig2 *tau2 *( J ))/{tau 1 * tau 1 *( 1 ) 
+tau2*tau2*(I)); 

eu=chi+et; 

pn->e=eu; 
pnp->e=et; 

alp=axl-2.*et; 

bet=bxl; 

gam=cxl-2.*et; 

pn->alp=(2.*bei+alp+gam)/2.; 

pn->bet=(bet-alp+gam); 

pn.>gain=(2.*bet-aIp+3.*gam)/2,; 

alp=ax2-2.*eu; 

bet=bx2; 

gam=cx2+2.*eu; 

pap->aip=(-2.*bet+aIp-gain)/2,; 
pnp->bet=alp+be t+gam ; 
pnp->gam=C2.*bet+alp+3.*gani)/2.; 

return (temp); 



double deriv(int a.int b.int c.int d,int p,int q)i 

double x,y,nep; 

x=p*p-q*q; 
y=-2*p*q; 
x=x/(p*p+q*q); 
y=y/(p*p+q*q); 

nep=:-(a*a+b*b+c*c+d*d)4-x*(a*a+b*b-c*c-d*d)-y*(2*a*c+2*b*d); 
retum(-27nep); 
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) 

struct complex bgener(struct edge *pc,struci edge *pn,double fc) 
( 

double fbar(int,ini); 
int a,b.c,d.bpd,apc; 
double e,aip,bet,gam; 
double axl.ayl,ax2,ay2; 
double bxl»byl,bx2»by2; 
double cx I ,cy 1 ,cx2»cy2; 
double sigl.iaul.sig2,tau2,chi; 
double vlc,vlu.vlt,v2c,v2u»v2t; 
double eu,et»fn.fnp; 
extern int count; 

double deriv(int,int,int,int.int.int); 
double dc.dn.dnp; 
struct complex temp; 
struct edge *pnp; 

pnp=pn; 
pnp++; 

a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 

apcs=a+c; 
bpd=lmi; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn-xl=a; 

pnp->a=d; 
pnp->b=-c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 

temp.x=rfn=fbar(a+apc,b+bpd) ; 

temp.y=fnp=:fbar(c+apc4+bpd); 

iemp.n=0; 
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couni+=^2; 

dc=deri v(d,-crb,a. 1.1); 
dn=deriv(d,-c,-b,a,2,l); 
dnp=deri v(d , - c b .a J ,2) ; 

fc=fc/dc; 
fn=fn/dn; 
fnp=fnp/dnp; 



axl=alp+4,*e; 
bxl=bet-4.*e; 
cxl=gam; 

ax2=alp+4.*e; 

bx2-bet; 

cx2=:gam-4,*e; 

chi=2.*(Cfc-gam-bet)/4.+e); 
) 



eise{ 

axl=alp+4*e; 

bxl=bet; 

cxl=gam+4*e; 

ax2=alp44*e; 
bx2=bet+4*e; 
cx2=gani; 

chi=2*((fc-gam-bet)/4-e); 
) 

V 1 c=-fn+cx 1 +(4. * bx 1 -3 . *ax 1 )/5- ; 

vlt=-4./5.; 

vlu=8./5.; 

v2c=-fnp+cx2+(4.*bx2+3.*ax2)/5.; 

v2u=4./5.; 

v2t=-875,; 

sigl=:vlc+chi*vlu; 
taul=:vlt+vlu; 

Sig2=v2c+chi*v2u; 
tau2=v2t+v2u; 



if(a>c){ 
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et=-(sig 1 *tau 1 *( 1 ) . _.g2*tau2*C 1 ))/(tau 1 *taLi I *(1 ) 
+tau2*tau2*(l)); 

eu=chi+et; 
pn->e=eu; 
pnp->e=et; 

alp=axl-2.*et; 

bet=bxl; 

gam=cxl-2,*et; 

pn->aip=(2 *bet+a!p+gam)/2.; 

pn->bei=(bet-alp+gam); 

pn->gam=(2,*bet-aIp+3.*gani)/2.; 

alp=ax2-2.*eu; 

bei=bx2; 

gam=cx2+2.*eu; 

pnp->aip-C-2.*bet+alp-gain)/2.; 

pnp->bet=alp+bet+gam; 

pnp->gam=(2.*bet+alp+3,*gam)/2.; 

return (temp); 
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//This is an implen...itation in C of a preferred embodiment 

//of the invention for the calculation of the inverse Fourier transform 



#include <stdio.h> 

#define SCAT 50 //SCAT=[exp sinh^(-J)(l/SCALE)]^2 

#define N 50 //input bandwidth 

#derine MINMODE 0 //minimum frequency of input 
#define PI 3.141592653589793 

#derine CHUNK 5000 //chunk of memory allocated as required 

struct complex { //a complex number 

double x; 
double y; 

); 

struct edge { //a complex arrow-structure 

int a,b,c,d; 

struct complex e,alp,bet,gam; 

}; 

struct pqxy{ //the basic data type of the output, where 

int p,q; //p,q corresponds to the complex number (p-iq)/(p+iq) 

double x.y; //at which the output function takes the complex value x+iy 
)» 

struct complex pfourier[2*N+l]; 

struct complex *fourier=pfourier+N; //the input Fourier coeffs 

//indexed -N -1,0,1....,N 

struct complex czer.cone.cmon ; //the modified 0.+ 1 1 Fourier coeffs 

struct pqxy graf[CHUNK]; //output spatial values in 

//clockwise-order around the circle 
//starting from the complex number -1 

//running tally of the total 
//number of output values 



//recursive routine for top of circle 
//recursive routine for bottom of circle 

//output normalization routine 
//input normalization routine 



int outcount=0; 



void main(void){ 



void genertop(struct edge *.int); 
void generbot(struct edge *.int); 

void normalout(void); 
void normalin(void); 



struct complex cogi(struct edge *); 



//returns complex wavelet coeff of argument 
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Struct complex ei,emu,emt,eu,ei: //wavelet coeffs on I,U^(-1 IT'^f-l )»U.T 

struct complex emtt,emtu,emut,emuu; //wavelet coeffs on T^{-2}, 

//r^{-l }U^{-J J.U'^f-l IT'^l-l },U^{-2) 

struct edge startfl]; //initial edge for each of 

//the various recursive calls 

int i; 

//begin input data 

for(i=0;i<=:2*N;i-H-)| 

pfourier[i].x=:0.: 
pfourier[i].y=0.; 
) 

fourier[5].x=.5; 
fourier[-5].x=.5; 

//end input data 



normalinO; //this stuffs czer,cone,cmon with the 

//modified 0,+l,-l Fourier coefficients 

graflouicount].p=:0; //initial normalized output value 

graf[outcount] .q= 1 ; 

graf[outcount] .x=0. ; 

graf[outcount] .y=0. ; 

outcount++; 



//calculate the wavelet coeffs for the various required edges of small generation 



start->a=start->d=l; 
start->c=0; 
start->b=-l; 
emt=cogi(start); 

start->a^tart->d=l ; 
start->c=0; 
start->b=l; 
et=cogi(start); 

start->a=start.>d=l; 
start->c=-l; 
start->b=0; 
emu=cogi(start); 

start->a=stan->d=l ; 
start->c=l ; 
start->b=0; 
eu=cogi(start); 

start->a=start->d= 1 ; 
start->b=start->c=0; 
ei=cogi(start); 

56 




wo 00/44104 



PCTAJS99/30584 



stan->a=start->d= 1 ; 
start->b=0; 
start->c=-2; 
emuu=cogi(start); 

start->a=start-xi= 1 ; 
start->b=-2; 
start->c=0; 
emtt=cogi(start); 

start->a=l; 
start->b=start->c=- 1 ; 
stan->d=2; 
emut=cogi(start); 

siart->a=2; 
start->d=l; 
start->b=siait->c=- 1 ; 
emtu=cogi(start); 



//serially perform the recursions for the initial edges in counter-clockwise order 



start.>a=start->d=l; //initial 

start->c=l; 

start->b=0; 

stan->e=eu; 

start->alp,x=2.*ei.x-2.*et.x; 

start->alp.y=2.*ei.y-2.*et.y; 

start->bet.x=2.*(emt.x-emu.x)-2.*ei.x; 

sian->bet.y=2.*(emt.y-emu.y)-2.*ei.y; 

siart->gam.x=2 *ei.x-2 *et.x; 

stan->gam.y=2.*ei.y-2.*et.y; 

genertopCstart, 1); //recursion for U 



start->c=0; 
start->b=l; 
siart->e=ei; 

start->alp.x=2*ei.x-2.*eu.x; 

start->alp.y=2.*ei.y-2.*eu.y; 

start->bet.x=2.*(ennt.x-emu.x)+2.*ei.x; 

start->bet.y=2.*(emt.y-eniu.y)+2.*ei.y; 

start->gam.x=-2.*ei.x+2.*eu.x; 

start->gam.y=-2.*ei.y+2.*eu.y; 

genertop(start,l); //recursion for T 



start->c=0; 

stan->b=-2; 

start->e=emtt; 

start->alp.x=:-(-2.*ei.x+2.*emu.x-4.*emt.x+2 *emut,x); 
start->alp.y=-(-2.*ei.y+2.*emu.y-4.*emt.y+2.*emul.y); 



stan->a=start->d=I; 



//initialize 



start->a=start->d=l; 



//initialize 
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start->bet.x=-(2.*c.....2.*emu,x+2.*emt.x); 
sian->bet.y=-(2.*ei,y-2.*emu.y+2.*emt.y); 
start->gam.x=2.*ei.x-2.*eniu.x+4.*emt.x-2.*emut.x; 
stan->gam.y=2.*ei.y-2.+emu.y+4.*emt.y-2.*emut.y; 



generbot(stari.l); 



//recursion forT^{-2) 



start->a=l; 
start->d=2; 



//initialize 



start->b=start->c=- 1 ; 
start->e=emut; 

start->alp.x=-(-2.*ei.x+2.*emu.x+2.*emt.x); 

start->alp.y=-(-2.*ei.y+2.*emu.y+2.*enit.y); 

siart->bet.x=-(2.*ei.x-2.*emu.x-6.*emt.x+4.*enm.x); 

start->bet.y=-(2.*ei.y-2.*emu.y-6.*emt.y+4.*emtt.y); 

start->gam.x=Z*ei.x-2.*emu.x-6*emt.x+4.*eintt.x; 

start->gam.y=2.*ei.y-2.*emu.y-6.*emt.y+4.+emtt.y; 

generbot(start,i); //recuision forU'^f-l )T^{-1 ) 



o.c«t //initialize 

start->d=l; 

siart->b=start->c=- 1 ; 

start->e=emtu; 

start->aIp.x=-(-2.*ei.x+2.*emu.x+2.*emt.x); 

start->alp.y=-(-2.*ei.y+2.*emu.y+2.*emt.y); 

siart->bet,x=-(-2.*ei.x+6.*emu.x+2.*enit.x-4.*emuu.x); 

start->bet.y=-(-2.*ei.y+6.*emu.y+2.*emt,y-4.*emuu.y); 

start->gam.x=-2.*ei.x+6.*emu.x+2.*emt,x-4.*emuu.x; 

start->gam.y=-2.*ei.y+6.*emu.y+2.*enit.y-4.*emuu.y; 

generbot(start,l); //recursion for T^{-1 )U^{-1 ) 



start->b=0; 
start->c=-2; 
siart->e=:emuu; 

start->alp.x=-(-2.*ei.x-4 *emu.x+2.*emt.x+2.*emtu.x); 

start->alp.y=-(-2,*ei.y.4.*eniu.y+2>emt.y+2.»emtu.y); 

start->bet.x=-(-2.*ei.x-2.*emu.x+2 *emt.x); 

start->bet.y=-(-2.*ei.y-2.*emu.y+2.*emt.y); 

start->gam.x=-2-*ei.x-4.*emu.x+2.*emt.x+2.*emtu.x; 

start->gam.y=-2.*ei.y-4.*emu.y+2.*emt.y+2.*emtu.y; 

generbot(start,l); //recursion for U'^{-2) 

//recursions complete, and graf is stuffed with normalized output data in counter-clockwise order 



start->a=stan->d= 1 ; 



//initialize 



normaloutO; 



//un-normalize by altering 
//each output value using the 
//modified 0.+l,-l Fourier coeffs 



printf("cmon=%e +(i) %e\n czer=%e +(i) %e\n cone=%e +(i) %e\n\n'\ 
cmon.x,cmon.y»czer.x,czer.y,cone.x,cone.y); 



for (i=0;i<outcount;i++) 
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printf("(%^. -^d) %e+(i)%e\n", 

grafTi].p,grar[i].q,graftiJ.x.graf[i].y); 



void genertop(struci edge *p.int g){ //recursion for top of the circle 

struct edge twofer[2],*ped; //recursion generates twofer[2] from p, where 

//the entries of twofer are the descendants of p 

struct complex temp; 

struct complex makecpx(int,int); //input p»q returns (p-iq)/(p+iq) 

double ct,st; //cos.sin for various angles 

void genertop(struct edge *»int); //recursive function 



void tgener(struct edge *,struct edge *); 



int subtendCstruct edge *).pp,qq; 



//tgener(p,twofer) stuffs twofer[2] 
//with the descendants of p 

//returns 1 if both argument edges 
//subtend angles less than SCALE, 
//otherwise returns 0 



tgener(p,twofer); 



if(subtend(twofer)) { 



ped=twofer; 
pp=-(ped->d); 
qq=^jed->c; 
grafIoutcount]-p=pp; 
grafloutcount].q=:qq; 
temp=makecpx(pp,qq); 
ct=temp.x; 
st=temp.y; 

graf[outcount].x=ped->alp.x 
+ped->gam.x-Kped- 

graffoutcount] .y=ped->alp.y 
+ped->gam.y+(ped- 

outcount++; 



//tgener(p,twofer) stuffs twofer 
//with the two descendants of p 

//if both descendants subtend 
//angles less than SCALE, then 
//write two output values and 
//terminate the recursion 

//stuff output values for the first 
//(counter-clockwise) sample point 



*c t+ped->bet.x * st 
•>e.x)*47(pp*pp+qq*qq); 
*'ct+ped->bet.y*st 
>e.y)*4./(pp*pp+qq*qq); 

//update number of output values 



//stuff output values for the second 
//(counter-clockwise) sample point 



ped-M-; 
pp=-(ped->d); 
qq=ped->c; 
graf[outcount] .p=pp; 
graf[ouicount].q=qq; 
temp=makecpx(pp,qq); 
ct=temp.x; 
st=temp.y; 

grafloutcount].x=ped->aIp.x*ct+ped->bet,x*st 

+ped->gam.x; 
graf[outcount].y=ped->aIp.y*ct+ped->bet,y*st 

+ped->gam.y; 

outcount++; //update number of output values 



return; 
) 



//terminate the recursion 



genertop(twofer,g+ 1 ); 



//in the contrary case that one of 
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genertop(twofer+l,£ . 
} 

void generboi(struct edge *p.int g){ 



struct edge twofer[2],*ped; 

struct complex temp, makecpx(int,int) 
double ct,st; 

void generbot(struct edge *,int); 

void bgener(struci edge *,struct edge *) 

int subtend(struct edge *),pptqq; 

bgener(p,twofer); 

if(subtend(twofer)) ( 



//the descendant edges subtends i. 
//large angle, then continue the 
//recursion in counter-clockwise sense 

//recursion for the bottom of 

//the circle, this is entirely 
//analogous to genertop. and 

//only the differences will 
//be noted here 



ped=:twofer, 

pp=-ped->b; //differs from genertop in that 

qq=ped->a; //b,a replaces d.c 

grafloutcount] .p=pp; 

graf[outcount].q=qq; 

temp=makecpx(pp.qq); 

ct=temp.x; 

st=temp.y; 

graf[outcount].x=-ped->alp.x*ct-ped->bet.x*st //differs from 

+ped->gam.x+(ped->e.x)*4./(pp*pp+qq*qq); // genertop in that 

grafIoutcount].y=-ped->alp.y*ct-ped->bet.y*st //alp and bet are 

+ped->gam.y+(ped->e,y)*4y(pp*pp+qq*qq); //multiplied by -1 

outcount++; 



ped-H-; 

pp=-ped->b; 

qq=I«d->a; 

graf[outcount].p=pp; 

grafloutcount] .q=qq; 

temp=makecpx(pp,qq); 

ct=temp.x; 

st=temp.y; 

grafIoutcount].x=-ped->alp,x*ct-ped->bei,x*st 

+ped->gam-x; 
graf[outcount].y=-ped->alp.y*ct-ped->bet.y*st 

+ped->gam.y; 
outcount++; 



//differs from genertop in 
//that b,a replaces d,c 



//differs from genertop 
//in that alp and bet are 
//multiplied by -1 



return; 
} 



generboi(twofer.g+ 1 ); 
generbot(i wofer-f 1 .g+ 1 ) ; 



struct complex mult{struct complex u.struct complex v){ //multiplication 

struct complex temp; //of complex numbers 
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temp,x=u.x*v,x-u.^ /.y; 
temp.y=u,x*v.y+v.x*u.y; 

retum(temp); 

} 



struct complex makecpx(int p,int q){ 
struct complex temp; 
double den; 

den=p*p+q*q; 

temp.x=(p*p-q*q)/den; 
temp.y=-2*p*q/den; 



) 



retuni(temp); 



int subtend(struct edge *p)( 
int nn; 

nn=(p->a)*(p->c)+(p.>b)*(p->d); 
if(nn*nn<SCAT) 
retum(0); 

P++; 

nn=(p->a)*(p->c)+(p->b)*(p->d); 

if(nn*nn<SCAT) 
retum(0); 

retimi(l); 



//input p.q returns 
//the complex number 
//(p-iq)/(p+iq) 



//subtend retums 1 if each of 
//of the two edges in the 
//array p subtends an angle 
//approximately less than 
//SCALE and retums 0 otherwise 



void normalin(void){ 

struct complex cz[4],co[4],cm[4]; 



struct complex temp; 
int n; 



cz[0].x=-l 
C2[0).y=0. 

C2[1].X=0. 

cz[l].y=0. 
cz[2].x=-l 
cz[2].y=0. 
cz[3].x=0. 
cz[3].y=0. 



co[0].x=0.; 
co[0].y=0.; 
co[l].x=-l, 
co[l].y=0.; 



//this routine calculates the modified 
//0,+l,-I Fourier coeffs czero,cone,cmon 
//from the input Fourier coeffs 



//intialize array for czero calculation 



//initialize array for cone calculation 
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cm[0].x=0.; 
cm(0].y=K).; 
cm[l].x=0.; 
cmi[l].y=0.; 
cin[2].x=0.; 
cm[2].y=-l.; 
cm[3].x=-l.; 
cm[3].y=0.; 



//initialize array forcmon calculation 



C2en=fourier[03; 
cone=fourier[l]; 
cmon=fourier[- 1 ]; 

for(n=2;n<=N;n++){ 

temp=mult(fourier[n],c2[n%4]); 

czer,x+=-temp.x; //contribution to czer from the 

czer.y+=-temp.y; , //positive Fourier coeffs 

temp=mult(fourier[-n],cz[(4*N-n)%4]); 



czer.x+=-temp.x; //contribution to czer from the 

czer.y+=-temp.y; //negative Fourier coeffs 

temp=:mult{fourier[n],co[n%4]); 

cone.x+=-temp.x; //contribution to cone from the 

cone.y+=-temp.y; //positive Fourier coeffs 

temp=muh(fourier[-n]»co[(4*N-n)%4]); 

cone.x+s-temp.x; //contribution to cone from the 

cone.y+=-temp.y; //negative Fourier coeffs 

temp=mult(fourierIn].cm[n%4]); 

cmon.x+=-temp.x; //contribution to cmon from the 

cmon.y+=-temp.y; //positive Fourier coeffs 

temp=mult(fourier[-n],cm[(4*N-n)%4]); 

cmon,x+=-temp.x; //contribution to cmon from the 

cmon.y+=-temp.y; //negative Fourier coeffs 



void normalout(void){ 



//this routine alters the 
//output values in the 
//array graf by adding 



ini n,m; 

struct complex iemp,templ,makecpx(int,int); 
double ct.st; 



//czero+cone*z+cmon*z'^(- 1 ), 
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//where z=(p-iq)/(pT . 

for (n=0;n<outcount;n-f+){ 

temp=makecpx(graf[n].p,graf[n].q); 
tempi =mult(temp,cone); 

temp.y=:-temp-y: 
temp=mult(temp.cmon); 

graf[n].x+=czer.x+templ.x+iemp.x; 
grafln] .y+=c2er.y+te mp 1 ,y +temp.y ; 

) 



struct complex cogi(struct edge *p){ //cogi(p) calculates the wavelet coeff for the 

//arrow underlying edge p from the Fourier coeffs 

int a»b»c.d»n; 

struct complex makecpx(int,int), mult(struct complex^struct complex); 
struct complex temp2,templ,temp,crun,cbl,cml,cb2,cm2; 

a=p->a; 

b=p->b; 

c=p->c; 

d=p->d; 

cnin.x=0.; 

crun.y=0.; 

cm 1 =makecpx(b,-a); 
cm2=makecpx(d,-c); 

cbl=CTil; 
cb2=cm2; 

for(n=2;n<=N;n-H-){ 

cbl=mult(cbl,cml); //updated nth power of cm! 
cb2=mult(cb2,cm2) ; //updated nth power of cm2 

if(n<MINMODE)//skip contribution to wavelet coeffs 

continue; //from low-frequency Fourier coeffs 

temp.x=n; 
temp.y=b*d+a*c; 
templ=mult(cbl,temp); 
temp,y=-temp.y; 
temp2=mult(cb2,lemp); 
temp.x=templ.x+temp2.x; 
temp.y=templ .y+temp2.y ; 
tempi=mult(fourier[n],temp); 
crun.x+stempl.x; 
crun.y+=templ.y; 

temp.x=-temp.x; 
temp2=mult(fourier[-n],temp); 

crun.x4-=temp2,x; //contribution to wavelet coeff from 

crun.y+=:temp2.y; //nth Fourier coeff for n negative 

) 



//contribution to wavelet coeff from 
//nth Fourier coeff for n positive 
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temp=crun; 

crun.x=-temp.y/4.; //overall multiple of i/4 

crun.y=temp.x/4.; 

retum(cnin); 

) 



void tgener(struct edge ♦pc,struct edge *pn)l 



//tgener(p,twofer) stuffs iwofer 

//with descendants of p for top 
//of circle 



int a.b,c,d,bpd.apc; 



//matrix entries a,b.c,d. 
/^pd=b-Ki and apc=a+c 



struct complex e,alp,bet,gam; //wavelet coefficient e, trig coefficients alp.bet,gam 



struct complex axl,ax2; 
struct complex bxl,bx2; 
struct complex cxl,cx2; 
int ayl,ay2,byl,by2,cyl,cy2; 

struct complex eu.et; 

struct complex temp; 



//index 1,2 refers to first,second 
//(counter-clockwise) descendants 
//prefix a,b.c refers to alp.bet,gam 
//update procedure below defines x,y 

//wavelet coeffs of first, second 
//descendant edges 



struct complex cogi(struct edge *),mult(struct complex^struci complex); 
struct complex makecpx(int.int); 

double ct.st; 

struct edge *pnp; 



pnp=pn; 
pnp-M-; 

a=pc->a; 
b=pc->b; 
c=pc->c; 
d=pc->d; 

apc=a4c; 
bpd=bKi; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 



//initialize pointers 
//initialize matrix entries 



//stuff matrix of first descendant 



//stuff matrix of descendant 
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alp=pc->alp; 
bet=pc->bet; 
gani=pc->gam; 



//initialize lagged trig coeffs 



e=pc->e; 



//initialize wavelet coeff 



if(a>c){ 



//prepare for stuffing alp.bet,gann in case a>c 



ax 1 .x=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; 
ax 1 .y=alp.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y ; 



ayl=2*(b*b-a*a); 

bxl.x=bet.x-Kc*d-a*d-b*c-a*b)*4*e.x; 
bxl.y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y; 

byl=4*a*b; 

cx 1 .x=gam.x+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.x; 
cx 1 .y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.y ; 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x-*4*(d*d-c*c)*e.x; 
ax2.y=alp.y-h4*(d*d-c*c)*e.y; 

ay2=2*(c*cKl*d); 

bx2.x=bet.x+8*e.x*c*d; 
bx2.y=bet.y+8*e.y*c*d; 

by2=-4*c*d; 

cx2.x=gam.x-4*e.x*(c*c-Kl*cl); 
cx2.y=gam.y-4*e.y*(c*c-Hl*d); 

cy2=2*(c*c4<i*d); 
) 

//prepare for stuffing aip.bet,gam in case c>=a 

ax 1 .x=alp.x44*e.x*(a*a-b*b); 
axl.y=:alp.y44*e.y*(a*a-b*b); 

ayl=2*(b*b-a*a); 

bx 1 .x=bet.x-8*e.x*a*b; 
bx I .y=bet-y-8*e.y *a*b; 

byl=4*a*b; 

cx 1 .x=gam.x+4*e.x*(a*a+b*b); 
cx 1 .y=gam.y+4*e.y *(a*a+b*b); 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; 
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ax2.y=ali.., r(a*a-b*t>f2*(b*ci-a*c)-c*c+d*d)*2*e.y; 
ay2=2*(c*c-d*d); 

bx2.x=bei.x+(a*ci-a*b+b*c+c*d)*4*e.x; 
bx2.y=bet,y+(a*d-a*b+b*c+c*d)*4*e.y; 

by2=-4*c*d; 



cx2.x=gam.x+(a*a+b*b- 
cx2.y=gam,y+(a*a+b*b- 

cy2=2*(c*c+d*d); 
} 

pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

pn->alp.x=ax 1 .x+ay 1 *et.x; 
pn->alp.y=axl.y+ayl*et.y; 

pn->bet.x=bx 1 -x+by 1 *et.x; 
pn->bet.y=bx 1 .y+by 1 ♦et.y; 

pn->gam-x=cx 1 .x+cy 1 ♦et.x; 
pn->gam,y =cx 1 .y+cy 1 ♦et.y; 

pnp->alp.x=ax2.x+ay2*eu.x; 
pnp->alp.y=ax2.y+ay2*eu.y; 

pnp->bet.x=bx2,x+by2*eu.x; 
pnp->bet.y=bx2.y+by2*eu.y; 

pnp->gam.x=cx2.x+cy2*eu.x; 
pnp->gam.y=cx2.y+cy2*eu.y; 



•2*(a*c+b*d)-c*c-d*d)*2*e.x: 
2*(a*c+b*d)-c*c-d*d)*2*e.y; 



//calculate wavelet coeffs 
//of the descendants 

//update first alp 
//update first bet 
//update first gam 
//update second alp 
//update second bet 
//update second gam 



} 



void bgener(struct edge *pc,struct edge *pn){ /^gene^(p.twofe^) stuffs twofer 

//with descendants of p for 
//bottom of circle, 

int a,b,c,d,bpd,apc; //it is entirely analogous to 

struct complex e.alp,bet,gam; //tgener, and only the 

//differences will be noted here 



struct complex axl,ax2; 
struct complex bxl,bx2; 
struct complex cxi,cx2; 

int ayl,ay2»byl»by2,cyl,cy2; 

struct complex eu,et; 
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Struct complex ten.^, 

struct complex cogi(struct edge *),mult(struct complex,struct complex); 
struct complex makecpx(int,int); 

struct edge *pnp; 

double ct.st; 



pnp=pn; 
pnp++; 

a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 



//differs from tgener in that d,-c»-b,a 
//replaces a,b,c,d 



apc=;a-K:; 
bpd=b+d; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 

pnp->a=d; 
pnp->b=-c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



//differs from tgener in that d,-c,-b»a 
//replaces a,b,c,d 



//differs from tgener in that drc,-b,a 
//replaces a,b,c.d 



if(a>c)[ 



axl.x=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; 
ax 1 .y=aip.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y ; 

ayl=2*(b*b-a*a); 

bx 1 .x=betx+(c*d-a*d-b*c-a*b)*4*e.x; 
bxl.y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y; 

byl=4*a*b; 

cx 1 .x=gam.x+(2*(a*c+b*d)-c*c-d*d+a *a+b*b)*2*e.x; 
cxl.y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.y; 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x-f4*(d*d-c*c)*e.x; 
ax2.y=alp.y-h4*(d*d-c*c)*e.y; 

ay2=2*(c*c-d*d); 

bx2.x=bet.x+8*e.x*c*d; 
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bx2.y=bei.^ .•8*e.y*c*d; 
by2=-4*c*d; 

cx2.x=gam.x-4*e.x*(c*c+d*d); 
cx2.y=gam.y-4*e.y*(c*c-Kl*d); 

cy2=2*(c*c+d*d); 
) 

else{ 

ax 1 .x=alp.x-»-4*e.x*(a*a-b*b); 
ax 1 .y=alp.y+4*e.y *{a*a-b*b); 

ayl=:2*Cb*b-a*a); 

bx 1 .x=bet.x-8*e.x*a*b; 
bxl .y=bet.y-8*e.y*a*b; 

byl=4*a*b; 

cxLx=gam.x44*e.x*(a*a+b*b); 
cxl.y=gam.y+4*e.y*(a*a+b*b); 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; 
ax2.y=alp.y+(a*a-b*b+2*(b*d-a*c>c*c+d*d)*2*e.y; 

ay2=2*(c*c-<l*d); 

bx2.x=bet.x+(a*d-a*b+b*c+c*d)*4*e.x; 
bx2.y=bety+(a*d-a*b+b*c+c*d)*4*e.y; 

by2=-4*c*d; 

cx2.x=:gam.x+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.x; 
cx2.y=:gam.y+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.y; 

cy2=2*(c*c+d*d); 
1 



pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

pn->alp.x=ax 1 .x+ay 1 *et.x ; 
pn->alp.y=ax 1 .y+ay 1 *et.y ; 

pn->bet.x=bx 1 .x+by 1 *et.x; 
pn->bet.y=bx 1 .y+by 1 *et.y; 

pn->gam.x=cx 1 .x+cy I *et.x; 
pn->gam.y=cx 1 .y+cy 1 *et.y; 

pnp->alp.x=ax2.x+ay2*eu,x; 
pnp->alp.y=ax2.y+ay2*eu.y; 

pnp->bet.x=bx2.x+by2*eu.x; 
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pnp->bet.y=bx2.y+u.. _*eu.y; 

pnp->gam.x=cx2.x+cy2*eu.x; 
pnp->gam,y=cx2.y+cy2*eu.y; 



//Here are the subroutines replacing tgener and bgener above in an 
//innplementaiion which takes advantage of renormalization 



void igener(struct edge *pc,struct edge *pn) 



int a,b,c»d,bpd,apc; 

struct complex e,alp,bei,gam; 

extern int outcount; 

extern struct pqxy graf[CHUNK]; 

struct complex axl,ax2; 

struct complex bxl,bx2; 

struct complex cxl.cx2; 

struct complex eu,et; 

struct complex cogi(struct edge *); 

double ct,st; 
struct edge *pnp; 

pnp=pn; 

pnp++; 

a=:pc->a; 

b=pc->b; 

c=pc->c; 

d=pc->d; 

apc=a+c; 
bpd=l>Ki; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 
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if(a>c){ 



ax 1 .x-alp.x+4.*e.x; 

ax 1 .y=alp.y+4.*e.y; 

bxl.x=bet.x-4.*e.x: 

bxl.y=bet,y-4.*e.y; 

cxl,x=gam.x; 

cxl,y=gam.y; 

ax2.x-axl.x; 

ax2.y=axl.y; 

bx2.x=bei.x; 

bx2.y=bet.y; 

cx2.x=gam.x-4.*e.x; 

cx2.y=gam.y-4.*e.y; 

} 



ax 1 .x=alp.x+4.*e.x; 

axl,y=alp.y+4.*e.y; 

bxl.x=bel.x; 

bxLy=bei.y; 

cx 1 .x=gam.x-f4. ♦e.x; 

cx 1 .y=gam.y+4.*e.y; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2,x=bet.x+4, *e,x ; 

bx2.y=bet.y +4. *e.y ; 

cx2.x=gain.x; 

cx2.y=gam.y; 

) 



pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

alp.x=axl.x-2.*et.x; 

bet.x=bxl.x; 

gain.x=cxl.x-2.*et.x; 

alp.y=axl.y-2.*et.y; 

betyssbxl.y; 

gani.y=cxl.y-2.*ei.y; 

pn->alp.x=(2.*bet.x+aIp.x+gani.x)/2.; 

pn->bet.x=(bet.x-alp.x+gam.x); 

pn->gam.x=(2.*bet,x-alp.x+3.*gam.x)/2.; 

pn->alp.y=:(2 *bet.y+aip.y+gam,y)/2.; 
pn->bet.y=(bet.y-alp.y+gam.y); 
pn->gam.y=(2.*bet.y-alp.y+3 ♦gam.y)/2.; 

alp.x=ax2.x-2.*eu.x; 

bet.x=bx2.x; 

gam.x=cx2.x+2-*eu.x: 

alp.y=ax2.y-2.*eu.y; 
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bei.y=bx2.y; 
gam.y=cx2.y+2.*eu.y; 

pnp->alp.x=(-2.*bei,x+aIp.x-gam-x)/2.; 

pnp->bet.x=alp.x+bet.x+gam.x; 

pnp->gam.x=(2.*bet.x+aIp,x+3.*gam.x)/2.; 

pnp->alp.y=(-2.*bet,y+aip.y-gam.y)/2.; 

pnp->bet.y=alp.y+bet.y+gam.y; 

pnp->gam.y=(2.*bet.y+alp.y+3.*gam.y)/2.; 

} 

double deriv(int ajnt b,int c,int d,int p.int q)( 

double x,y,nep; 

x=p*pHq*q; 
y=-2*p*q; 
x=x/(p*i>+<)*q); 
y=y/(p*p+q*q); 

nep=-(a*a+b*b+c*c-Ki*d)+x*(a*a+b*b-c*c-d*d)-y*2.*(a*c+b*d); 
reium(-2-/nep); 

) 

void bgener(struct edge *pc,struct edge *pn) 
i 



int a.b.c,d,bpd,apc; 

struct complex e,alp,bet,gain; 

extern int outcount; 

extern stmct pqxy graflCHUNK]; 

struct complex axl,ax2; 

struct complex bxl»bx2; 

struct complex cxl,cx2; 

struct complex eu,ei; 

struct complex cogi(struct edge *); 

double ct,st; 

struct edge *pnp; 

pnp=pn; 
pnp-M-; 



a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 

apc=^+c; 
bpd=b-Ki; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 
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pnp->a=d; 
pnp->b=-c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 



if(a>c){ 

ax 1 .x=alp.x+4.*e.x; 
axl.y=alp.y+4,*e.y; 
bx 1 .x=bet.x-4.*e.x; 
bx 1 .y=bet.y-4.*e.y ; 
cxLx=gam.x; 
cxl.y=gam.y; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2,x=bet.x; 

bx2.y=bet.y; 

cx2.x=gam.x-4.*e.x; 

cx2.y=gam,y-4.*e.y; 

) 



else{ 

axl.x=alp.x+4.*e.x; 

ax 1 ,y=alp.y+4.*e.y; 

bxl.x=bet.x; 

bxl.y=bet.y; 

cx 1 .x=gam.x+4. *e.x; 

cxl.y=gam.y+4.*e.y; 

ax2.x=:axLx; 

ax2.y=axl.y; 

bx2.x=bet.x+4.*e.x; 

bx2.y=bet.y+4.*e.y; 

cx2.x=gam.x; 

cx2.y=gam.y; 

) 

pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

pn->e=eu; 
pnp->e=et; 

alp.x=axl.x-2.*et.x; 

bet.x=bxl.x; 

gam.x=cxl.x-2.*et.x; 

alp.y=axl.y-2.*ei.y; 
bet.y=bxl.y; 
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gam.y=cxl.y-2.*t. ^ , 

pn->alp.x=(2.*bet.x+alp.x+gam.x)/2.; 

pn->bet.x=(bet.x-alp.x+gam.x); 

pn->gam.x=(2.*bet.x-alp.x+3.*gani.x)/2.: 

pn->alp,y=(2.*bet.y+aIp,y+gam.y)/2.; 

pn->bet.y={bet.y-alp.y+gam.y); 

pn->gam.y=:(2.*bet.y-aip.y+3.*gam.y)/2.; 

alp.x=ax2.x-2.*eu.x; 
beLx=bx2.x; 
gam.x=cx2.x+2. *eu.x; 

alp.y=ax2.y-2.*eu.y; 

bet.y=bx2.y; 

gam.y=cx2.y+2,*eu.y; 

pnp->a]p.x=(-2.*bet.x+alp.x-gam.x)/2.; 

pnp->bet.x=alp.x+bet.x+gam.x; 

pnp->gam.x=(2.*bet.x+aIp.x+3.*gam.x)/2.; 

pnp->alp.y=(-Z*bet.y4-alp.y-gam.y)/2.; 

pnp->bet.y=alp.y+bet.y+gam.y; 

pnp->gam.y=C2.*bet.y+alp.y+3.*gam.y)/2,; 
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Methods of D al Filtering and Multi-Dimensional Date ompression 
Using the Farey Quadrature and Arithmetic. Fan, and Modular Wavelets 

Claims 

I claim: 

1. A method of compressing and reconstructing input data comprising the steps of: 
providing input data in an initial form; sampling the input data at the points of the 
Farey quadrature; transforming the Farey-sampled data into an intermediate form; 
quantizing the Farey-sampled data while in the intermediate form; using the quan- 
tized intermediate form to store or transmit the Farey-sampled data; de-quantizing 
the quantized intermediate form of the Farey-sampled data prior to reconstructing the 
data; and reconstructing the intermediate form of the Farey-sampled data in order to 
obtain reconstructed data for use. 

2. A method as recited in claim 1 wherein: the step of transforming the Farey-sampled 
data into an intermediate form is achieved via a wavelet transformation; and the inter- 
mediate form of the Farey-sampled data takes the form of wavelet coefficients. 

3. A method as recited in claim 2 wherein the wavelet transformation is based on 
arithmetic wavelets. 

4. A method as recited in claim 3 wherein the arithmetic wavelets are defined in 
accordance with the following expressions: 



MO) = 



-f 2 cos ^ + 2 sin 5 - 2 
+2 cos ^ - 2 sin ^ + 2 
-2 cos ^ - 2 sin 6-2 
-2 cos 5 + 2 sin e + 2 



if 0 <d < f, 

if f < 5 < TT, 

iin <0 
if ^ < ^ < 27r; 



and if^=^° d^'^^O l)"''^ labels an arithmetic arrow; 

and if A labels a top arithmetic arrow and a > c, then 

' +4(^2 _ c2)cos e + 8cdsine- 4(rf2 + c^); 

if tan-^gl^ < 9 < tan-^ (,y,ff(ffj, , 

-f-2(a2 + d2 _ ^2 _ ^2 ^ 2ac - 2bd) cos 6 
+4(cd — ab-bc- ad) sin 9 
+ 2(a2 + 62 _ c2 _ rf2 + 2ac + 2bd); 

if tan-T<2ig^g±^ < ^ < tan- V 



■dAie) = { 



26a 



+2(62 + d2 - a2 - c2 + 2ac - 2bd) cos 9 
+4{ab + cd-ad- be) sin 6 
+ 2(-a2 _ fe2 _ ^2 _ ^2 ^ 2ac + 2bd); 

if tan->j;|^ <9< tan-^.?^f-''>l"-°l 



(d-by-'-(c-a)'' ' 



1^0; 



otherwise; 
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and if A labels a top i hmetic arrow and c > a, then 



( +2(a2 -h c2 - 62 _ ^2 _ + 2bd) cos 6 
+4(ad + fee — a6 - cii) sin 5 
+2(a2 + 62 _^ ^2 _^ ^2 _ 2ac - 26d); 

It tan < 6^ < tan '^yf^, 

+2(a2 4- rf2 „ 52 _ ^2 ^ 26d - 2ac) cos 6 
+4(ad -hbc-hcd- ab) sin 6 
-f 2(a2 + 62 _ c2 - d2 - 26d - 2ac); 

II tan 53r^ S fc' S tan it-^^y^l\a-^c)^^ 

if tan"^ 2(b4-/)(a+c) < ^ ^ fon-l 

»^ 0; otherwise; 
and if A labels a bottom arithmetic arrow and a > — c, then 

( +2(a2 + c2 - 62 - 4- 2ac - 26rf) cos 6 
H-4(-ad -be- ab- cd) sin ^ 
+2(a2 4- 62 -I- c2 + d2 ^ 2ac + 26d); 

if 4. — -1 2{b+d){a+c) ^ n ^ 4. -1 26a 

« tan ^^^p^j^y-^-—^. < 6^ < tan '55^, 

+2(62 + c2 - a2 - d2 4- 2ac - 2bd) cos 0 
4-4(a6 - ad - 6c - cd) sin ^ 
4- 2(c2 4- d2 - a2 - 62 4- 2ac + 26d); 

iftan-^^<0<tan-^^|i|^^^ 

+4 (c2 - d2) COS e - 8cd sin 9 + 4(c2 4- d^); 

if tan~^ 2(a--bKc-a) <^ ^ ^ ^.^--i 2cd 
" (d-6)2-(c-a)i ^ €^ S tan 

^ 0; otherwise; 
and if A labels a bottom arithmetic arrow and — c > a, then 

^ 4-4(62 - a2) cos ^ -h Sab sin 5 - 4(a2 + 62); 

if tan-^^ < e < tan-^^^I^p^i^ 

+2(62 + c2 - a2 - d2 + 26d - 2ac) cos 0 
+4(ad + 6c + (26 - cd) sin 5 
+ 2(c2 + d2 - a2 - 62 - 2ac - 26d); 

^a{0) = ^ if ^^^-'ii^W^^ < G < tan-^^, 

+2(62 + d2 - a2 - c2 + 26d - 2ac) cos 6 
+4 (ad + 6c + cd + a6) sin 6 
+ 2(-a2 - 62 - c2 - d2 - 2ac - 26d); 

if tnn~^ 2crf < ^ < 4- — -1 2(64-rf)(a-4-c) 
"tan ^TT^T S S tan -^f^^^^^.^^^^^i , 



^9^(5) = < 



U; 



otherwise. 
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5. A method as reci in claim 2 wherein the wavelet transfc .ation is based on 
modular wavelets. 

6. A method as recited in claim 5 wherein the modular wavelets are defined in accor- 
dance with the following expressions: 



n>0 



n>0 



where 



^T-» [9) + 2 sin 6; if = T-^, 



V).^=(J O'-'^^'O O'"^ 



+2 COS ^ + 2 sin ^ - 2 

•d (9^ = I "'■^ cos e - 2 sin 5-1-2 

' 1 -2 cos 5 - 2 sin 5 - 2 

-2 cos 5 -f 2 sin 5 -t- 2 



if 0 < 5 < f , 

if f < e < TT, 

if TT < e < ^, 
if ^ < 5 < 27r; 



and if^— d^'^^O l)""^ labels an arithmetric arrow; 
and if A labels a top arithmetic arrow and a > c, then 



^a{6) = < 



' +4(^2 _ c2)cos ^ -I- 8cd sin e - 4.{dP + c^); 

if tan-ig^ < e < t^n- \^i';/l%l%, , 

+2(a2 + -b^-c'^ + 2ac - 2bd) cos 6 
-H4(cd — a6 — be - ad) sin 5 
+ 2{o? + 62 - c2 - d2 ^. 2ac 4- 26£i); 

if tan-^ (,2^yH°X:), < e < tan-^5^, 

-f 2(62 -f- ^2 _ ^2 _ ^2 ^ 2oc - 26d) cos 6 
+A{ab + cd-ad- be) sin 9 
+ 2(-a2 - 62 - c2 - d2 2ac -I- 26d); 

if tan-i^ < e < tan- \^[iyI'}_%-_% , 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 
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' +2(a2 + c2 - &2 _ ^2 _ + 2bd) cos 5 
+A{ad + bc — ab~ cd) sin 5 
+2(a2 -f- 62 ^ ^2 ^ ^2 _ _ 26d); 

if t^n~^ 2(d-b)(e-a) ^ a ^ 2c£i 
" (d-6)i-(c-ap < f < tan ^533^, 

+2(a2 + d2 - 62 _ c2 + 26<i - 2ac) cos 5 
,9 /fl-v _ } +4(ad + 6c + cd - a6) sin 6 

^^)-\ + 2(a2 + f,2 _ ^2 _ ^2 _ 26d _ 2ac); 

if tan-3lS<e<tan-i^gM5^i^, 

+4(a2 - 62) cos e - Sab sin 6 + 4(a2 + 62); 

" '^^'^ (6+d)i!-(o+c)!2 < f < tan 
^ 0; otherwise; 
and if A labels a bottom arithmetic arrow emd a > — c, then 



' +2(a2 + c2 - 62 - rf2 + 2ac - 26d) cos 6 
+4(—ad ~bc — ab — cd) sin 6 
+2(o2 + 62 + c2 + d2 ^ 2ac 4- 26rf); 

;f * — -1 2(ft+d)(o+c) ^ n ^ . -1 26o 

" ^^"^ (6+d)Ma+c)2 <6'<tan 'e^tst, 

+2(62 + c2 - a2 - d2 4- 2ac - 26d) cos 6 
+-4(a6 — ad — be — cd) sin d 
+ 2(c2 + d2 - a2 - 62 + 2ac + 26d); 

if tan-^5^^ < 0 < tan-^ JiyHpJ^^ , 

+4 (c2 - d2) cos ^ _ gj„ Q ^ ^ ^2). 

" tan < < tan 'grr^, 

' 0; otherwise; 
and if A labels a bottom arithmetic arrow amd — c > a, then 

' +4(62 - a2) cos 5 + 8a6 sin d - 4(a2 + 62); 

if tan~^ ^2&a < < <■ — -1 2(d-6)(e-a) 

+2(62 + c2 - a2 - d2 + 26d - 2ac) cos 0 
+4(ad + 6c + a6 — cd) sin 9 
+ 2(c2 + d2 - a2 - 62 - 2qc - 26d); 
^a{6) = <! if tan- \2£y4°4, < 0 < tan-^^^, 

+2(62 + d2 - o2 - c2 + 2bd - 2oc) cos 0 
+4(ad + 6c + cd + a6) sin 0 
+ 2(-a2 - 62 - c2 - d2 - 2ac - 26d); 

if t.in~' 2cd < <- «. — -1 2(b+d)(o+c) 



10; 



otherwise. 
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7. A method as recite .1 claim 2 wherein the wavelet transform. :>n is based on fan 
wavelets. 

8. A method as recited in claim 7 wherein the fan wavelets are defined in accordance 
with the following expressions: 



n>0 



n>0 



where 



i{A^U-\T-\ 



^T-^iO) + 2 sin 9; i{A = T-\ 



with S = 



r +2 cos e + 2 sin 0-2 

+2 cos ^ - 2 sin 6 + 2 

-2 cos ^ - 2 sin e - 2 

L -2 cos 61 + 2 sin ^ + 2 



if 0 < e < f , 
if f < ^ < TT, 
if TT < ^ < ^, 
if ^ < 0 < 27r; 



andifyl=^° d}'^(o l)"''^ labels an arithmetric arrow; 
and if A labels a top arithmetic arrow and a > c, then 



' +4(^2 _ c2)cos 5 + 8cd sin 0 - 4(d^ + c^); 

if tan-^5^ < 0 < tan-^^|Ig;^K^, 

+2(o2 +cP-b^-c^ + 2ac - 2bd) cos 6 
+4(cd - a6 - 6c - ad) sin 6 
+ 2{a? + 62 _ c2 _ rf2 ^ 2ac + 26<i); 

+2(62 + d2 - a2 - c2 + 2ac - 26d) cos 9 
+4(a6 + cd-ad- be) sin ^ 
+ 2(-a2 - 62 - c2 - d2 + 2ac + 26d); 

if tan-'^a^ < 9 < tan-^ JiyH— 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 
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' +2{a? + 0"^ -b"^ -cf - 2ac + 2bd) cos 6 
+4(ad + bc- ah- cd) sin 6 
+2{o? + 62 + c2 + d2 _ 2ac - 2bd); 

+2(o2 + d2 - 62 - c2 + 26d - 2ac) cos 9 
+4(ad + 6c + cd - a6) sin 6 
+ 2(a2 + 62 - c2 - d2 - 26d - 2ac); 



2ed 



if ta.n~^ 2ed < /Q < t^^-l 2(fc+d)(a+e) 

+4(a2 - 62) cos 9 - Sab sin 9 + 4(a2 + 62); 

V 0; otherwise; 
and if A labels a bottom arithmetic arrow and a > — c, then 

^ +2(a2 + c2 - 62 _ ^2 ^ 2ac - 26rf) cos 0 
+4(-ad - 6c - a6 - erf) sin 0 
+2(a2 + 62 + c2 + + 2ac + 26d); 

if tan~^ 2(6+d)(a+c) <^ ^ ^ ton-l 26a 
" (6+d)ii-(a-+-c)2 ^ ^ S tan p3^) 

-f 2(62 _^ ^2 - a2 - rf2 ^ 2ac - 2bd) cos 5 
,9 (a\ _ ^ +4(a6 — ad — 6c — cd) sin ^ 

- "I + 2(c2 + d2 - a2 - 62 + 2ac + 26d); 

if tan->^ < 9 < t^n- \^^_^,^l%-_% , 
+4 (c2 - d?) COS ^ - 8cd sin e + 4(c2 + ^2); 



0; otherwise; 
and if A labels a bottom arithmetic arrow and — c > a, then 



' +4(62 - a2) cos ^ 4- 8ab sin 6 - 4(^2 -|- 62); 

II tan ^^3^ s S tan ^— p-ii—j^, 

+2(62 + c2 - a2 - ^2 ^ 2bd - 2ac) cos 0 
+4 (ad + 6c -f a6 - cd) sin 0 
+ 2(c2 + d2 - a2 - 62 - 2ac - 26d); 



^a(0) = < 



" (d-6)3l(c-a)2 ^ ^ S tan ^s_ca , 

+2(62 + d2 _ a2 - c2 + 26d - 2ac) cos 0 
+4(ad + 6c + cd + a6) sin 0 
+ 2(-a2 _ 62 - c2 - d2 - 2ac - 26d); 

It tan ^T3^<^<tan ^^^j^^^-^-^ , 



lO; 



otherwise. 
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J- A method as recitec claim 2 wherein the steps of sampling tl nput data at the 
points of the Farey quadrature and transforming the Farey-sampled data into wavelet 
coefficients is characterized by a binary cascade of arrow structures arising from the 
calculation of wavelet coefficients by sampling the input data in its initial form at the 
points of the Farey quadrature as defined along a circle in the complex plane. 

10. A method as recited in claim 2 wherein the wavelet transformation is based on 
arithmetic wavelets and the step of transforming the Farey-sampled data into arithmetic 
wavelet coefficients includes the step of regularization. 

11. A method as recited in claim 2 wherein the wavelet transformation and the recon- 
struction are based on a finite number of wavelets. 

12. A method as recited in claim 2 wherein: wavelet transformation is achieved via a 
plurality of wavelet transformation operations each for a specified arc along a circle in 
the complex plane; and the plurahty of wavelet transformation operations are calculated 
in parallel with each other. 

13. A method as recited in claim 1 wherein the quantized intermediate form of the 
data is stored and retrieved prior to de-quantizing and reconstructing the data for use. 

14. A method as recited in claim 1 wherein: the steps of transforming the Farey- 
sampled data into the intermediate form and quantizing the intermediate form of the 
Farey-sampled data occur at a first physical location; the quantized intermediate form 
of the data is transmitted to a second physical location; and the steps of de-quantizing 
the quantized intermediate form of the Farey-sampled data and reconstructing the 
intermediate form of the Farey-sampled data in order to obtain reconstructed data 
occurs at the second physical location. 

15. A method as recited in claim 14 further comprising the step of transmitting a 
reconstruction algorithm to the second physical location in order to accomplish recon- 
struction at the second physical location. 

16. A method as recited in claim 1 wherein the input data is digital data. 

17. A method as recited in claim 1 wherein the input data is analog data and the 
analog input data is sampled by placing analog sensors at locations corresponding to 
the points of the Farey quadrature. 

18. A method as recited in claim 1 wherein the input data and the reconstructed data 
are multi-dimensional. 

19. A method as recited in claim 2 wherein the step of reconstructing the intermediate 
form of the Farey-sampled data in order to obtain the reconstructed data is achieved 
via an inverse wavelet transformation which is characterized by a binary cascade. 

20. A method as recited in claim 2 wherein the calculation of the wavelet transformation 
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includes the step of rei malization. 

21. A method as recited in claim 1 wherein the sampling of data at the points of the 
Farey quadrature includes temporary storage of sequences of data in a buffer. 

22. A digital signal processing method for obtaining data in a transform domain com- 
prising the steps of: providing input data in an initial form; sampling the input data at 
the points of the Farey quadrature; using an intermediate transform to transform the 
Farey-sampled data into an intermediate form; and using a primary transform to con- 
vert the data in intermediate form of the Farey-sampled data into transform coefficients 
representative of the data in a transform domain. 

23. A method as recited in claim 22 wherein the primary transformation consists of 
one of the groups of following transformations: Fourier, Hilbert, Haar, Laplace, Bessel, 
Laguerre, Hermite, Chebyshev, Hotelling, Mersenne and Fermat. 

24. A method as recited in claim 22 wherein the primary transformation consists of the 
Fourier transform, and the intermediate transform is a wavelet transformation based 
on arithmetic wavelets. 

25. A method as recited in claim 24 wherein the intermediate transform is an arithmetic 
wavelet transform and the primary transform is a Fourier transform, and the conversion 
of the data from the intermediate form to the Fourier coefficients, c^, is given by the 
following expressions: 





and the coefficients Cq^, , cf j are given by 



(a-c)2zp22(6-d)(a-c)]/2 
{a + c)2 =F 2z(6 -f d)(a -f- c)]/2, 



^4 = 29h{c''-{-£ 



) - 9r.[{b-df^{a^cf] 



where the angles are 



Oh = tan~^( 



2c6 



), e^ = tan-i( 



2(6 -frf)(a -he) 



(b + d)2 - (a + c)2 
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26. A method as recited in claim 22 wherein the intermediate transform is a wavelet 
transform based on arithmetic wavelets, and the calculation of the intermediate trans- 
form includes the step of renormalization. 

27. A method as recited in claim 22 wherein the primary transformation consists of the 
Fourier transform, and the intermediate transform is a wavelet transformation based 
on modular wavelets. 

28. A method as recited in claim 22 wherein the intermediate transform is a wavelet 
transform based on modular wavelets, and the calculation of the intermediate transform 
includes the step of renormalization, 

29. A method as recited in claim 22 wherein the primary transformation consists of the 
Fourier transform, and the intermediate transform is a wavelet transformation based 
on fan wavelets. 

30. A method as recited in claim 22 wherein the intermediate transform is a wavelet 
transform based on fan wavelets, and the calculation of the intermediate transform 
includes the step of renormalization. 

31. A method as recited in claim 22 wherein the sampHng of data at the points of the 
Farey quadrature includes temporary storage of sequences of data in a buffer. 

32. A method as recited in claim 22 wherein the input data is multi-dimensional. 

33. A method as recited in claim 22 wherein the constants and coefficients required 
in the calculation of the intermediate transform and in the calculation of the primary 
transform are archived and recovered from memory at run time. 

34. A method of data processing comprising the steps of: providing input data in an 
initial form; sampling the input data at the points of the Farey quadrature; transform- 
ing the Farey-sampled data to a second form; and processing the data in the second 
form. 

35. A method as recited in claim 34 wherein: the step of transforming the Farey- 
sampled data into the second form is achieved via a wavelet transformation; and the 
second form of the Farey-sampled data takes the form of wavelet coefficients. 

36. A method as recited in claim 35 wherein the wavelet transformation is based on 
arithmetic wavelets. 

37. A method as recited in claim 36 wherein the arithmetic wavelets are defined in 
accordance with the following expressions: 
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4-2 cos e -i- 2 sin (9 - 2; if 0 < ^ < |, 

9 {fi\ _ I +2 cos ^ - 2 sin 0 + 2; if f < 5 < tt, 

^/ W - S _2 cos ^ - 2 sin 5 - 2; if tt < g < 

-2 cos 0 -I- 2 sin e + 2; if ^ < 0 < 2;?; 



and if^=^^ ^^^^^ labels an arithmetric 

and if A labels a top arithmetic arrow and a > then 



arrow: 



( +4(d2 _ c2)cos e + 8cd sin 9 - A{d? + c^); 

II tan gs^r^ » ^ tan (6^rf)2_];a+c)S > 

+2(a2 + - 62 _ c2 4- 2ac - 2bd) cos 6» 
+4(cd — at — 6c — ad) sin ^ 
+ 2{a? + b'^-c^ -d^ + 2ac + 26d); 

if t.in-^ 2(fe+d)(a+c) < « <. 26a 

+2(62 + rf2 - a2 - c2 + 2ac - 26d) cos & 
+4(a6 + cd — ad — be) sin d 
+ 2(-a2 _ 62 _ c2 - d2 4. 2ac + 26d); 

if tan-ij^ < ^ < tan-^ Jiynpj^^ , 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 



f +2(a2 + c2 - 62 _ ^2 _ 2ac + 26d) cos 6 
-\-4{ad + 6c — a6 — cd) sin ^ 
-{-2(a2 + 62 -h -f d2 _ 2ac - 26rf); 

if fnn-^ 2(d-b)(c-a) ^ p ^ j. -1 2cd 

" (d-f,)2-(c-a)2 S < tan 573;-?, 

+2(a2 4. rf2 _ ^2 _ ^2 ^ 26£f - 2ac) cos 9 
-\-4{ad + 6c -h cd — a6) sin 6 
+ 2(a2 4- ft2 _ ^2 „ ^2 _ 2bd - 2ac); 

" tan s S tan (^^^^^^(a+c)^ > 

-h4(a2 - 62) cos e - Sab sin 5 -f- 4{a^ + 62); 

0; otherwise; 



and if A labels a bottom arithmetic arrow and a > —c, then 
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' +2(a2 + c- -b--d'^ + 2ac - 2hd) cos 6 
+4(-ad -bc-ab- cd) sin d 
+2(a2 + 62 c2 + d2 + 2ac + 2bd); 

if tan-^ 2(fe+d)(a+c) < ^ fan-l 26a 

+2(62 + c2 - o2 _ d2 ^. 2ac - 26d) cos 5 
+4(a6 - ad-hc — cd) sin 6 
+ 2(c2 + ^2 _ o2 _ 62 ^ 2ac + 26d); 

+4 (c2 - d2^ COS e - 8cd sin ^ + 4(c2 + d^); 

"tan (<i_t)ii_\c-a)i < f < tan ^gsr^r, 



^ 0; otherwise; 
and if labels a bottom aritiimetic arrow and — c > a, then 



' +4(62 _ q2) ^,q3 ^ g^j gjjj ^ _ 4^^2 + ^2). 

" tan bj_„s i: w 5- tan (d-fc)^T^e-tt)a > 

+2(62 + c2 - a2 - d2 ^. 2bd - 2ac) cos 6 
+4(arf + 6c + a6 - cd) sin ^ 
+ 2(c2 + ^2 - a2 _ 62 - 2ac - 26d); 

if tan-^ Jlt7Mc-S'^ ^ ^ ^ tan-i^^, 



+2(62 + d2 - o2 - c2 + 26d - 2ac) cos 6 
+4(ad + 6c + cd + ab) sin 5 
+ 2(-a2 _ 62 _ p2 _ ^2 _ - 26d); 

if tan-'^^f^ <e< tan-i7#f^i;^iiS±4 



(6+d)^-(a+c)*' 



otherwise. 



38. A method as recited in claim 35 wherein the wavelet transformation is based on 
modular wavelets. 



39. A method as recited in claim 38 wherein the modular wavelets are defined in 
accordance with the following expressions: 



n>0 



n>0 

where 

f i?A(e); \iA^u-\T-\ 

[ i?7'-'(<9) +2 sin 61; if = r-\ 
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with S 



r +2 COS 0 + 2 sin 0 - 2 

+2 cos e - 2 sin 0 + 2 

-2 cos 5 - 2 sin <9 - 2 

I -2 cos 5 + 2 sin <9 + 2 



ana 



if 0 < 0 < f , 

if f < (9 < TT, 

if TT < ^ < ^ , 
if ^ < 0 < 27r; 



and if ^=^^ ^^^^^ l)"^"^ ^^^^^^ arithmetric 
and if >1 labels a top arithmetic arrow and a > c, then 



arrow; 



^ +A{dP - c2)cos 6^ + 8cd sin 6* - 4(^2 -f- c^); 

iftan-^^<0<tan-^^fg-^ 

+2(a2 + rf2 _ _ ^2 _^ 2ac - 26d) cos 5 
4-4 (cd -ab-bc" ad) sin 5 
4- 2(a2 + 62 _ ^2 _ ^2 ^ 2ac + 2bd); 

^AiO) = ^ tan-^^^^^^ < e < tan"^^, 

+2(6^ + rf2 - a2 „ ^2 _^ 2ac - 26d) cos 6 
-f-4(a6 H- cd - ad - 6c) sin 6 
+ 2(-a2 - 62 - c2 - d2 + 2ac -h 26d); 

II tan 5!n-r S c/ S tan ^5„b)5T(c-a)5 > 

^ 0; otherwise; 
and if A labels a top arithmetic arrow and c > a, then 



' +2(a2 + c2 - 6^ - d2 _ + 26d) cos 0 
+4(ad + 6c - a6 - cd) sin ^ 
-f 2(a2 4- 62 + c2 + rf2 - 2ac - 26d); 

if ^r-^-^ 2(d-b)(c-~a) ^ /) ^ 2cd 

+2(a2 + - 62 - c2 + 26d - 2ac) cos 6 
4- 4(ad 4- ic 4- cd — a6) sin 0 
4- 2(a2 4- 62 » c2 - d2 - 26d - 2ac); 

if tan-i^ < e < t^n-^^^i^^^S^ 

4-4(a2 _ 62) cos 0 - 8ab sin 0 4- 4(o2 4- b^); 

if tan~^ 2(6+d)(a+c) <^ ^ ^ tnn~l 26a 



otherwise; 



and if A labels a bottom arithmetic arrow and a > — c, then 
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( +2(a2 + c- -b'^ -d- + 2ac - 2hd) cos 6 
+A{-ad -be - ab- cd) sin 9 
+2{a? + 62 + c2 + d2 + 2ac + 2bd); 

if t.in-^ 2(b+d)(a+c) ^ a ^ ton-l 26o 

" '^^^ (6+d)i'-(a+c)a < f < tan ^Erf^T, 

+2(62 + c2 - _ ^. - 26d) cos 6 
+4(o6 -ad -be- cd) sin 5 
+ 2(c2 + d2 - a2 - 62 + 2ac + 26d); 

if tan-^ < 9 < tan-^ ^,y)(--g, . 

+4 (c2 - d2) cos 5 - 8cd sin 5 + 4(c2 + cP); 

" '^an (d26)^-(c-a)ii < e < tan ^^f^^, 

0; otherwise; 
and if A labels a bottom arithmetic arrow and — c > a, then 

" ^an s c S tan ^^-jp^iL—t^^ 

+2(62 + c2 - a2 - d2 + 26d - 2ac) cos 6 
+4(ad + 6c + a6 - erf) sin ^ 
+ 2(c2 + d2 _ ^2 _ ^2 _ 2ac - 26d); 

^a{Q) = { if tan-^ ^.^lffH— < e < tan-15^, 

+2(62 + d2 - o2 - c2 + 26d - 2ac) cos 9 
+4(ad + 6c + cd + 06) sin 6 
+ 2(-a2 - 62 - c2 - d2 _ 2ac - 26d); 

if tan-^^ < e < tan-^ .^"i^l^li. , 



otherwise. 



40. A method as recited in claim 35 wherein the wavelet transformation is based on 
fan wavelets. 



41. A method as recited in claim 40 wherein the fan wavelets are defined in accordance 
with the following expressions: 



n>0 



n>0 

where 

■>-i 



^a{9) = < i9u-i{e)-2sine; if A = U-\ 
l^T-^(e)+2sme; ifA = T-\ 
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with S 



I andT= (J J );ana 



f +2 cos (9 + 2 sin 6> - 2; if 0 < ^ < f , 

+2 cos ^ - 2 sin e + 2; if f < 6* < tt, 

-2 cos 5 - 2 sin 6^ - 2; H tt < 0 < 

I -2 cos /9 -f 2 sin ^ + 2; if ^ < 0 < 27r; 



and if^=^^ ^^^^^ l)^^ labels an arithmetric 
and if A labels a top arithmetic arrow and a > c, then 



arrow; 



' 4-4(^2 - c2)cos 0 + Serf sin ^ - 4(d^ + c^); 

if tan"^-^;^ < 9 < tan-^ -2(^i^)i£±£)_ 
n can ^^r:-? S S tan ^^^j-^^^-^-—^ , 

+2(a2 + ^2 _ ^2 _ ^2 ^ 2ac - 26(i) cos 0 
+4(ccf -ab-bc- ad) sin 5 
-f 2(a2 4- 62 „ ^2 „ ^2 _^ 2ac + 26d); 

= i if tan-^^^igi^i^ < ^ < tan"^^, 

+2(62 + - a2 „ ^2 ^ 2ac - 2bd) cos 0 
+4(ab -hcd — ad — be) sin 5 
+ 2(~a2 - 62 - c2 - rf2 ^ 2ac + 26d); 

if tan-^^ < ^ < t.n-^^M^^(^^ 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 



' +2(a2 + - 62 - rf2 _ -h 2bd) cos (9 
+4 (ad + 6c - a6 — cd) sin 5 
+2(a2 4- 62 + c2 + d2 _ 2ac - 26d); 

if ^^^''xi^^^^ < e < tan-^^, 

+2(a2 + d2 _ j,2 _ ^2 ^ 26d - 2ac) cos 0 
i9 — ) +4(ad -f- 6c H- cd — a6) sin 0 

- ^ + 2(a2 4- 62 - c2 - d2 - 26d - 2ac); 

if tan~^ 2cd < n i. -1 2(6+d)(a^-c) 



+4(a2 - 62) cos 0 - 8a6 sin <9 + 4(a2 + 62); 
w 0; otherwise; 



if tan"^ . 2(6+d)(a+c) < ^ ^-^-,-1 26a 



and if A labels a bottom arithmetic arrow and a > — c, then 
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( +2(a2 -f c2 - 62 _ + 2ac - 2bd) cos 6 
+4{—ad -be — ab- cd) sin 9 
+2(a2 + ^2 4. c2 + + 2ac + 26d); 

:f f nn-^ 2(fr+d)(a+e) ^ a ^ -1 2ba 

" (6+d)-i-(<i+c)2 S f < tan Err^, 

+2(62 ^. c2 _ a2 _ ^2 ^ 2ac - 26<i) cos 6 
+4(a6 - ad — 6c — cd) sin 6 
+ 2(c2 + d2 _ ^2 _ ^2 ^. + 26d); 

if tan-i^ < e < tan-^ (,^<yHpaW 

+4 (c2 - d2) COS ^ - Scd sin ^ + 4(c2 + d2); 

if tan-^ < e < tan-^5|^, 

k 0; otherwise; 
and if A labels a bottom arithmetic arrow and — c > a, then 

' +4(62 _ (;os 9 + 8a6 sin 6 - 4(a2 + 62); 

if tan-i^^ < e < tan-^ J(t^/)(— 

+2(62 + c2 - a2 - d2 4. _ 2ac) cos (9 
+4(ad + 6c + a6 - cd) sin ^ 
+ 2(c2 + d2 _ a2 _ 62 _ - 26d); 

= { < e < tan-i^lf^, 

+2(62 + d2 _ q2 _ ^2 ^ 26d - 2ac) cos 9 
+4(ad + 6c + cd + 06) sin 6 
+ 2(-a2 - 62 - c2 - d2 _ 2ac - 26d); 

if tan-^^l^ < 9 < tan-^ ^yKffg, , 



Otherwise. 



42. A method as recited in claim 35 wherein the steps of sampling the input data at 
the points of the Farey quadrature and transforming the Farey-sampled into wavelet 
coefficients is characterized by a binary cascade of arrow structures arising from the 
calculation of wavelet coefficients by sampling the input data in its initial form at the 
points of the Farey quadrature as defined along a circle in the complex plane. 

43. A method as recited in claim 42 wherein the wavelet transformation is based on 
arithmetic wavelets and the step of transforming the Farey-sampled data into arithmetic 
wavelet coefficients includes the step of regularization. 

44. A method as recited in claim 35 wherein the wavelet transformation is based on a 
finite number of wavelets. 

45. A method as recited in claim 35 wherein: wavelet transformation is achieved via a 
plurality of wavelet transformation operations each for a specified arc along a circle in 
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the complex plane; anc e plurality of wavelet transformation opei ons are calculated 
in parallel with each other. 

46. A method as recited in claim 34 wherein the input data is digital data. 

47. A method as recited in claim 34 wherein the input data is analog date and the 
analog input data is sampled by placing analog sensors at locations corresponding to 
the points of the Farey quadrature. 

48. A method as recited in claim 34 wherein the input data is multi-dimensional. 

49. A method as recited in claim 34 wherein the sampling of data at the points of the 
Farey quadrature includes temporary storage of sequences of data in a buffer. 

50. A digital processing system comprising a digital filter that receives an input signal 
and outputs a transformed signal, the digital filter being an arithmetic wavelet filter in 
which the arithmetic wavelets are defined in accordance with the following expressions: 



+2 cos 0 + 2 sin 0 - 2 

4-2 cos 5 - 2 sin 5 -h 2 

-2 cos 5 - 2 sin 6-2 

-2 cos e + 2 sin 0 4- 2 



if 0 < e < f , 

if f < 6> < TT, 

if 7r<(9< 

if ^ < 61 < 27r; 



and ifA=^^ ji: l)'^''^ labels an arithmetric 
and if A labels a top arithmetic arrow and a > c, then 



arrow; 



' +4(^2 - c2)cos e + 8cd sin 0 - 4{<f H- c^); 

if tan-^^ < 0 < tan-^ ^yifc^g. , 

+2(a2 4- cJ2 - 62 _ c2 + 2ac - 2bd) cos 0 
-f 4(cd — ab-bc — ad) sin 0 
4- 2(a2 4- 62 _ c2 - rf2 ^ 2ac + 2bd); 

d^{0) = { if tan-^^lgi^i^ < 0 < tan-^fe, 

4-2(62 + d2 - a2 - c2 + 2ac - 2bd) cos 0 
4-4(a6 -hcd — ad — be) sin 0 
+ 2(-a2 - 62 - c2 - d2 4- 2ac 4- 26d); 

if tan-^^ < 0 < tan-^^I|^^^^ 



lO; 



otherwise: 



and if A labels a top arithmetic arrow and c > a, then 
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C +2{o? + c'^~b'^-d- - lac + 2hd) cos 6 
+4 (ad -Vhc- ah — cd) sin 0 
+2(a2 + fe2 + ^2 + ^2 _ - 2M); 

+2(a2 + d2 _ ^2 _ ^2 + 26d - 2ac) cos 6 
+4(ad + bc + cd- ab) sin d 
+ 2(a2 + 62 - c2 - d2 - 2bd - 2ac); 

if tan-^3|^ < 9 < tan-^^^gi^i^. 

+4(a2 - 62) cos 9 - Sab sin 9 + 4(0= + 62). 

if * — -1 2(6+d)(tt+e) ^ a ^ 26a 

"tan < C < tan 'j;^, 

>> 0; otherwise; 
and if ^ labels a bottom arithmetic arrow and a > — c, then 



' 4-2(a2 + c2 - 62 _ ^2 ^. - 2bd) cos 6 
+4(-ad - 6c - 06 - cd) sin ^ 
+2(a2 + 62 + c2 + rf2 + 2ac + 2bd); 

+2(62 4- c2 - o2 - d2 + 2oc - 2bd) cos ^ 
+4(a6 — ad — 6c — cd) sin ^ 
+ 2(c2 + d2 - a2 - 62 + 2ac + 26d); 

if tan-^^ < 6 < tan-^ Jiy?{— 4, , 

+4 (c2 - d2) cos 5 - 8cd sin 9 + 4(c2 + d2). 

if t^n- \^%^-J^^, < 9 < tan-ig^s^, 

>. 0; otherwise; 
and if A labels a bottom arithmetic arrow and — c > a, then 

+4(62 _ ^2) cos 9 + 8a6 sin 9 - A{a^ + 62); 

if tan-^s^ < e < tan- \2(y2i— J). . 

+2(62 + c2 - a2 - d2 + 26d - 2ac) cos 9 
+A{ad + bc + ab — cd) sin 9 
+ 2(c2 + d2 - a2 - 62 - 2ac - 26d); 
^a{0) = I if tan- \,^(t7^^-), < ^ < tan-^s^, 

+2(62 + d2 _ ^2 _ c2 + 26d - 2ac) cos 9 
+4 (ad + 6c + cd + a6) sin 9 
+ 2(-a2 - 62 - c2 - d2 - 2oc - 26d); 

if tan-15^ < 9 < t^n-^ ^,Z^.fl%X%, , 



lO; 



otherwise. 
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51. A digital processir ystem comprising a digital filter that rec 3s an input signal 
and outputs a transformed signal, the digital filter being a modular wavelet filter in 
which the modular wavelets are defined in accordance with the following expressions: 



n>0 



n>0 



where 



^a{6) = < du-iie)- 2 sin 6; i{A=^U-\ 
[dT-iiO)+ 2 sin 9; ifyl = T-^ 



r -1-2 cos ^ + 2 sin e - 2 

+2 cos 0 - 2 sin 0 + 2 

-2 cos 5 - 2 sin e - 2 

I -2 cos ^ + 2 sin 5 -f 2 



if 0 < <9 < f , 

if f < e < TT, 
if TT <^< ^, 

if ^ < ^ < 27r; 



and if>l=^^ ^^^^^ l)^"^ ^^^^^^ arithmetric 
and if A labels a top arithmetic arrow and a> c, then 



arrow; 



-h4(<P - c2)cos 0 + 8cd sin 0 - 4{d'^ + c^); 

iftan-^^<^<tan-^^i|^^ 

-h2(a2 -i- d2 _ _ ^2 ^ 2ac - 2bd) cos 0 
+4(cd — ab — be — ad) sin 0 
+ 2(a2 + f,2 _ ^2 _ ^2 ^ 2ac -h 26d); 

if trtn-^ 2(64-d)(a+c) n ^ ^ -1 26a 

" tan (6+d)2„(a4-c)2 < < tan ptr^, 

-h2(62 + d2 „ ^2 _ ^2 _^ 2ac - 2bd) cos 6 
+4(a6 -\- cd — ad — be) sin 5 
-h 2(-a2 ^ 62 _ c2 - d2 ^ 2ac + 26rf); 

" tan ^73^ s (7 s tan (d,fc)-iT(c-a)a * 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 
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( +2(a2 + c2 - 62 _ ^2 _ 2ac + 2bd) cos 6 
+4{ad + bc — ab — cd) sin 6 
+2(a2 + 62 ^ ^2 + rf2 _ 2ac - 2bd)\ 

:f 2(d-6)(c-a) a ^ 2ed 

+2(a2 + d2 - 62 - c2 + 26d - 2ac) cos ^ 
9 /a\ — J +4(ad + bc+ cd — ab) sin 6 
^A(t^) ~ <j ^ 2(a2 + 62 _ c2 - d2 _ 26d - 2ac); 

ir tan 373^ s f S tan (6+^)^,(0^^)'' ' 
+4(a2 - 62) cos 5 - Sab sin ^ + 4(a2 + 62); 

k 0; otherwise; 
and if A labels a bottom arithmetic arrow and a > — c, then 



( +2(o2 + c2 - 62 - d2 + 2ac - 26d) cos 6 
+4(—ad — be — ab — cd) sin 6 
+2(o2 + 62 + c2 + d2 + 2ac + 26d); 

if t.iTi-^ 2(6+d)(a+c) -a 4. -1 2bo 

" (6H-d)'-i-(a+c)S S f S tan prrjT, 

4-2(62 + c2 - a2 - d2 + 2ac - 26d) cos 6 
+4(a6 — ad — bc — cd) sin 5 
+ 2(c2 + d2 _ ^2 _ ^2 ^. 2ac + 2bd); 

if tan-^s^ < e < tan-^ Jiy)<— 

+4 (c2 - d2) COS ^ - 8cd sin ^ + 4(c2 + ^2). 

" *a.n (j_j)a_);c-a)4 ^ ^ tan grr^r, 

k 0; othei*wise; 
and if A labels a bottom arithmetic arrow and — c > a, then 

( +4(62 _ q2) cos ^ + 8ab sin 5 - 4(a2 + 62); 

if tan-^5^^ < e < tan- \,2(y2<— 

+2(62 + c2 _ a2 _ ^ ^. 26d _ 2ac) cos ^ 
+4 (ad + 6c + a6 — cd) sin 0 
+ 2(c2 + rf2 _ ^2 _ 62 _ 2ac - 26d); 

dm = ^ if tan-^ ^21ff-(:-a)^ ^ ^ ^ tan-i^lE^, 

+2(62 + d2 _ ^2 _ c2 + 26d - 2ac) cos ^ 
+4(od + 6c + cd + o6) sin 6 
+ 2(-a2 _ 62 - c2 - d2 _ 2ac - 26d); 

if tan-^5|^<e<tan-^^|i|i^^g±g,, 



Otherwise. 
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62. A digital processin, /stem comprising a digital filter that recc ;s an input signal 
and outputs a transformed signal, the digital filter being a fan wavelet filter in which 
the fan wavelets are defined in accordance with the following expressions: 



n>0 



n>0 



where 



i9a W ={du-^ {6) - 2 sin ^; ilA = U-\ 



»uh5=(; j),a„dr=(j iy..n. 



J- +2 cos ^ 2 sin ^ - 2; if 0 < ^ < f , 

+2 cos 5 - 2 sin ^ + 2; if f < 0 < tt, 

-2 cos <9 - 2 sin 6> - 2; if tt < (9 < 

I -2 cos e + 2 sin 6^ + 2; if ^ < 5 < 27r; 



and if^— ^^^^^ labels an arithmetric arrow; 

and if A labels a top arithmetic arrow and a > c, then 



' -f-4(d2 - c2)cos e -f- Serf sin 0 - 4(d^ + c^); 

iftan-^5l!^<0<tan-^^^i|5^^ 

+2(a2 + d2 _ 52 _ ^2 ^ 2ac - 26d) cos 6 
+4(cd — ab—bc — ad) sin ^ 
+ 2(a2 + 62 _ ^2 _ ^2 ^_ 2ac -f 2M); 

= { tan-^^^gi^i^ < ^ < tan-^^, 

-h2(62 + rf2 _ ^2 _ ^2 _^ 2ac - 2bd) cos 6* 
+4(a6 + cd - ad - 6c) sin ^ 
+ 2(-a2 - 62 „ ^2 _ ^2 _^ 2ac + 26d); 

1^ ^an s tan 



otherwise; 



and if A labels a top arithmetic arrow and c > a, then 
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+2{a^ + c- -b"^ - 2ac + 2bd) cos 6 
+4 (ad + bc-ab- cd) sin 6 
+2(a2 + fe2 ^ c2 + d2 - 2ac - 2bd); 

if tnn~^ 2(rf-b)(c-a) ^ n ^ 4. -1 2erf 

+2(a2 + d2 _ ^2 _ ^2 ^ 2bd - 2ac) cos 6 
+4(ad + 6c + cd - 06) sin ^ 
+ 2(a2 + 62 _ g2 _ ^2 _ _ 2ac); 

if tan-^^l^ < ^ < tan-i^^ig^I^, 
+4(a2 - 62) cos d - 8a6 sin 6 + 4(0^ + 62); 

if <-nn-l 2(6+d)(a+e) .^^ /i ^ , -i 26a 

- 0; otherwise; 
and if A labels a bottom arithmetic arrow and a > —c, then 



' +2(a2 + c2 - 6^ - d2 ^ 2ac - 2bd) cos 6 
+4(-ad - 6c - a6 - cd) sin ^ 
+2(a2 + 62 + ^2 ^ ^2 ^ 2ac + 26d); 

if tan- \/^S'-^i:::>. < ^ < tan-^^^, 

+2(62 + c2 - a2 - d2 ^ 2ac - 26d) cos 9 
,9 //?\ _ J +4(a6 — ad — be — cd) sin ^ 
^'^i W - <j ^ 2(c2 + d2 - a2 - 62 + 2ac + 26d); 

if tan-^sl^ < 5 < tan-i 
+4 (c2 - d2) cos e - Scd sin 6 + 4{c^ + (P); 

if tan~^ 2(d-b)(c-a) < < tnn-l-2cd 

w 0; otherwise; 
and if ^ labels a bottom arithmetic arrow and — c > o, then 



' +4(62 _ q2) cos d + 8a6 sin 6 - A{a^ + 62); 

if tan-i^ < e < t^n- \^[iy!'l[Z% , 

+2(62 ^ c2 - a2 - d2 ^ 26d - 2ac) cos 0 
+4(ad + 6c + a6 - cd) sin 6 
+ 2(c2 + d2 _ ^2 _ ^2 _ - 2bd); 

+2(62 + c{2 _ q2 _ c2 + 2fed _ 2ac) cos 9 
+4 (ad + 6c + cd + a6) sin 9 
+ 2(-a2 _ 62 - c2 - d2 - 2ac - 26d); 

if tan-'J^ <9< tan-^ , !'ittf^''V\, ■ 



otherwise. 
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53, A method of obtai .g an inverse transform of input digital a comprising the 
steps of: providing input digital data in an initial form; transforming the data in the 
initial form to an intermediate form; and transforming the intermediate form of the 
data using an intermediate inverse transform to obtain output values at the points of 
the Farey quadrature. 



54. A method as recited in claim 53 wherein the input digital data is in the form 
of Fourier coefficients and the combination of transforming the data from the Fourier 
coefficients to the intermediate form and from the intermediate form to obtain output 
values at the points of the Farey quadrature constitutes an inverse Fourier transform. 

55. A method as recited in claim 54 wherein the intermediate form of the data is in 
the form of wavelet coefficients. 



56. A method as recited in claim 55 wherein the wavelet coefficients are coefficients for 
arithmetic wavelets, 

(^a(0); if yl^f/-^T-^ 

^A(d) = < ?9t/-i(e) + 2sin&; ilA = U-\ 
[ T^^-i (0) - 2 sin e; if >1 = T-i, 

with C/~^ ^ (-1 l}' "^"^ (o 1^)' conversion of the data from the 
initial form to the intermediate form is based on the following expressions: 



Arte 



= cos 710 -\- i sin nO 



where the sum is over all arithmetic arrows, the underlying chord of which has complex 
endpoints ^, 77, and 



Cq — 



-1, 


n 


= 0(4) 


0, 


n 


= 1(4) 


-1, 


71 


= 2(4) 


0, 


71 


= 3(4) 



( 0, n = 0(4) 

-1, n = l(4) 

i, n = 2(4) 

I 0, n = 3(4) 



0, n = 0(4) 

0, n = 1(4) 

-i, n = 2(4) 

L -1, n = 3(4). 



57. A method as recited in claim 55 wherein the wavelet coefficients are coefficients 
for modular wavelets, and conversion of the data from the initial form to the 
intermediate form is based on the following expressions: 



e*"* = cos ne + i sin nO 



= -[c^ + c^e'^ + cHie-'*] 

+ i +T?") + ^ {C - ^")|[V'J/^(^) - 2V>i(5) + ^t;-M(^)], 
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where the sum is over <• irithmetic arrows, the underlying chord c hich has complex 
endpoints ^,77, where U^^ (±1 j)' ^"^^ "^^^ere 



-1, n = 0(4) 

„ ^ I 0, n = 1(4) 

° ] -1, n = 2(4) 

0, n = 3(4) 



Cj — 



0, n = 0(4) 
-1, n = l(4) 

1, n = 2(4) 
0, n = 3(4) 



0, n = 0(4) 

0, n = 1(4) 

-i, n = 2(4) 

I -1, n = 3(4). 



58. A method as recited in claim 55 wherein the wavelet coefficients are coefficients for 
fan wavelets, and conversion of the data from the initial form to the intermediate 
form is based on the following expressions: 

e'"* = cos nO + i sin ne 

+ i E{"(r + + ^ (C - ^")} [</>a(0) - 

where the sum is over all arithmetic arrows, the underlying chord of which has complex 
endpoints ^, tj, where ^ = (^l , and where 

— 

Ci — 





n 


= 0(4); 




n 


= 1(4); 




n 


= 2(4); 


I 0, 


n 


= 3(4); 



0, 


n = 0(4); 




r 0' 


n = 0(4); 


-1, 


n = 1(4); 




0, 


n = 1(4); 


i, 


n = 2(4); 






n = 2(4); 


0, 


n = 3(4); 




l-l, 


n = 3(4). 



59. A method as recited in claim 53 wherein the intermediate form of the data is in 
the form of wavelet coefficients. 

60. A method as recited in claim 53 wherein the initial form of the input digital 
data is Fourier coefficients, the intermediate form of the data is in the form of wavelet 
coefficients, and the output values at the points of the Farey quadrature are obtained 
via the use of a binary cascade of arrow structures arising from the calculation of 
wavelet coefficients as defined along a circle in the complex plane. 

61. A method as recited in claim 60 further comprising the step of storing terminal 
arrow structures in order for restart or iterative refinement. 

62. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for 
arithmetic wavelets. 

63. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for 
modular wavelets. 

64. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for 
fan wavelets. 



65. A method as recited in claim 59 wherein the input data is multi-dimensional. 
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FIGURE la 

(b+d)+i(a+c) 



(b+d)-i(a+c) 




FIGURE lb 



c 




(d-b)+i(c-a) 



(d-b)-i(c-a) 

FIGURE Ic 




FIGURE 2 
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(b-h2d)+i(a+2G7 
(b+2d)-i(a+2c) d+ic 





FIGURE 1 1 
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Start 



Routine: main 

Data, Variables, and Effects: 

input: data presented as a function f(p,q), 
giving the function value at the 
point (p-iq)/(p-i-iq) 
output: array of complex Fourier coefficients 
effect: calculates Fourier transform 



1. Normalize input: 

a) Initialize arrray of complex fourier coefficients to zero 

b) Compute abar, bbar, cbar so that 

-2 pq 



2 2 

— p -q 

f(p,q)=f(p,q)- [abar* — : + bbar* 



p2+q2 



p^+q2 



-fcbar] 



is zero for (p,q)=(l,0), (0.1), and (1.1) 
Then: f(p,q) returns normalized values above 



I 



2. Start recursive calculation of arrows on top of 
the circle with call to gener (Figure 6) with initial 
values: arrows doe. generation=l,envy=0, f(-l,l) 



I 



3. Start recursive calculation of arrows on bottom of 
the circle with call to gener (Figure 6) with initial 
values: arrows doe, generations 1, en vy=0, f(l,l)=0 



4. Normalize output: divide Fourier coefficients by pi 



5> Display computed coefficients and exit 



FIGURE 5 
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15. End cascade and make 
correction calling prune 
(Figure 8) with arrow as 
argument 



± 



recursion 
x2 ^ 




7. Generate two descendant 
arrows in cascade using 
least-squares regularization 



8. Update Fourier coefficients 
calling Charles (Figure 7) 
with each descendant arrow 



14. Recursively call gener 
with the two descendant 
arrows, envy <-- envy, 
generation<--generation+ 1 



Routine: gener 

Data, Variables, and Effects: 

input: arrow, generation, envy, function value 
output: void 

effect: recursive routine to generate arrows and 
update Fourier coefficients 




FIGURE 6 
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Start 



1 6. Clear flag 



I 



17. Calculate preliminary 
quantitites 



Routine: Charles 

Data. Variables, and Effects: 

input: one arrow 

output: flag representing whether array 

update for arrow is non-negligible 

effect: modifies array of Fourier coefficients 

in accordance with the formulas in Step 2 
for the Fourier transform 




19. Pass data to routine for 0, 
+1,-1 Fourier coefficients 




27. Proceed to next 
Fourier coefficient 





24. Compute contribution of 
current arrow to current 
Fourier coefficient 



FIGURE 7 
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Start 

T 



28. Calculate final updates taul,tau2 of the two lagged 
trigonometric functions for terminal processing as at 
the end of the discussion of Fourier transform 



29. Calculate trigonometric extrapolations sigmal,sigma2 
for terminal processing as at the end of the discussion 
of Fourier transform using 3 further data samples for each of 
the two descendant arrows 



Routine: prune 

Data, Variables, and Effects: 

input: one arrow 
output: void 

effect: modifies array of Fourier coefficients 

by implementing the fmal processing of 
arrows in the cascade 



30. Modify the Fourier coefficients by adding those of 
sigmal-taul and sigma2-tau2 truncated as at the end 
of the discussion of Fourier transform 




32. Pass data to routine for 0, 
+1,-1 Fourier coefficients 



33. Return 



X 
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34. Calculate modified Fourier coefficients 



m 



c„ 
n m 



ngtO,-l,+l 



form=0,-l,+ l, where the coefficients c 
are given at the beginning of the section 
on inverse Fourier transform 



m 



35. Using the formulas from Step 1 at 
the beginning of the section on inverse 
Fourier transform, calculate the wavelet 

coefficients for A = I, U, T, U "\ T * ^, 

y-lj-l .p-ly-l 



36. For A=U. T, T U"^ T^^ 



use the wavelet coefficients from program 
segment 2 to stuff an arrow labeled A 
with the correct lagged trigonometric 
function and wavelet coefficient and 
serially call gener (Figure 10) with these 
arrows as argument thus starting these 
six recursive calculations of arrows 



I 



37. Correct the output function values by 
adding + c^j e"^"^ + c^^ ei^S 



I 



38. Display computed function values 
and exit 



Routine: main 

Data, Variables, Effects: 

input: array of complex Fourier coefficients 
c for Inl <= N 

output: array of three-tuples (p,q,z), where 
output function takes complex value 
2 at (p-iq)/(p-f-iq) given in the counter- 
clockwise order on the circle 

effect: calculates inverse Fourier transform 



FIGURE 9 
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Start 



recursion 
x2 



A 



41. Recursively call gener with the 
two descendant arrows as argument 



Routine: gener 

Data, Variables, and Effects: 

input: one arrow 
output: void 

effect: recursive routine to generate 
arrows and compute function 
values 



39. Generate the two descendant arrows 
in the cascade calculating arithmetic 
wavelet coefficients using the formulas 
from Step 1 from the section on inverse 
Fourier transform 




42. Perform the final update of 
lagged trigonometric functions 



43. Stuff output array with function values 
at endpoints of the two descendant arrows 



44. Return 







FIGURE 10 
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